In $\mathbb{R}^n$ we have some coordinates, so let $x^i (t)$ be the coordinate presentation of the curve. A shorthand notation $\dot{x}^i (t_0) = \frac{d}{dt}x^i |_{t=t_0}$ simplifies the calculations, as usual. The geodesic equation then says that $ (\nabla_{\dot{x}}\dot{x})^k = \frac{d^2 x^k}{dt^2} + \dot{x}^i \dot{x}^j \Gamma^k_{ij} = 0 \tag{1} $ where $x(t)$ has a unit speed parametrization.
Recall that the Christoffel symbols are given by $ \Gamma^k_{ij} = \frac{1}{2} g^{kl} (g_{il,j} + g_{lj,i} - g_{ij,l}) \tag{2} $
Remark. In the standard Euclidean metric $g_{ij} = \delta_{ij}$ in $\mathbb{R}^n$ the Christoffel symbols vanish and equation (1) gives straight lines as its solutions.
Let's now look what happens with the Christoffel symbols if we replace $g$ with $\hat{g} = e^{2 \rho} g$. Observe that $ \hat{g}_{ij,l} = 2 \rho_l e^{2 \rho} g_{ij} + e^{2 \rho} g_{ij,l} $ and substitute this into (2) to get $ \begin{align} \hat{\Gamma}^k_{ij} &= \frac{1}{2} e^{-2 \rho} g^{kl} ( 2 \rho_j e^{2 \rho} g_{il} + e^{2 \rho} g_{il,j} + 2 \rho_i e^{2 \rho} g_{lj} + e^{2 \rho} g_{lj,i} - 2 \rho_l e^{2 \rho} g_{ij} - e^{2 \rho} g_{ij,l}) \\ &= \Gamma^k_{ij} + \rho_j \delta^k{}_i + \rho_i \delta^k{}_j - \rho^k g_{ij} \tag{3} \end{align} $ where we have used metric $g$ to raise index $k$ in $\rho_k := \partial_k \rho$
This transformation affects equation (1) in the following way: $ \begin{align} (\hat{\nabla}_{\dot{x}}\dot{x})^k &= \frac{d^2 x^k}{dt^2} + \dot{x}^i \dot{x}^j \hat{\Gamma}^k_{ij} \\ &= \frac{d^2 x^k}{dt^2} + + \dot{x}^i \dot{x}^j (\Gamma^k_{ij} + \rho_j \delta^k{}_i + \rho_i \delta^k{}_j - \rho^k g_{ij}) = 0 \end{align} $
Your equation is now obtained immediately from the above calculations since you start with the Euclidean metric $g_{ij} = \delta_{ij}$ for which the Christoffels vanish, as has been already noted, so you only need to take into account that $\operatorname{grad}(\rho)\cdot \dot{x} = \rho_i \dot{x}^i$ and $\dot{x} \cdot \dot{x} = \dot{x}^i \dot{x}_i$
Ah, well, $\operatorname{grad}(\rho)^i = \rho^i = g^{ij} \rho_j$, of course.
Indeed, the above considerations are very formal and straightforward. In order to get a deeper insight one needs to contemplate the conformal transformations of $\mathbb{R}^n$. In dimension $n=2$ the situation is controlled by the complex analysis, while in dimensions 3 and higher the Liouville's theorem states that all such transformations are compositions of translations, dilations, and inversions (so called Moebius transformations).
(The Einstein summation convention is used throughout this post)