I don't know what \rho_0 might be.
The Lagrangian for the electromagnetic field (minimally) coupled to the (necessarily conserved) electromagnetic current reads (using Heaviside-Lorentz units with c=1 and the west-coast convention for the metric, \eta_{\mu \nu}=\text{diag}(1,-1,-1,-1)
\mathcal{L}=-\frac{1}{4} F_{\mu \nu} F^{\mu \nu}+J_{\mu} A^{\mu}.
Now you can go one step further and write down the Lagrangian for a charged particle and the electromagnetic current coupled to the electromagnetic field.
The free-particle part reads
L=-m \sqrt{1-\dot{\vec{y}}^2},
and the current density is given by
j^{\mu}(x)=q \int_{\mathbb{R}} \mathrm{d} t \; \dot{\vec{y}} \delta^{(4)}(x-y).
Here, \vec{y}(t) is the trajectory of the particle as a function of the coordinate time in the inertial reference frame with y^0=t.
Note that the total action is a Lorentz scalar and thus the equations of motion are relativistically consistent although not written in manifest covariant way. The action now consists of three parts:
The free action of the em. field
S_{\text{f}0}[A]=\int \mathrm{d}^4 x \left (-\frac{1}{4} F_{\mu \nu} F^{\mu \nu} \right ),
the free action of the particle
S_{\text{p0}}[\vec{y}]=-m \int \mathrm{d} t \sqrt{1-\dot{\vec{y}}^2},
and the interaction term
S_{\text{int}}=\int \mathrm{d}^4 x A_{\mu} j^{\mu}=\int \mathrm{d} t A_{\mu}[t,\vec{y}(t)] \frac{\mathrm{d} x^{\mu}}{\mathrm{d} t}=\int \mathrm{d} t \left [q A^0(t,\vec{y}(t))-\frac{\mathrm{d} \vec{y}}{\mathrm{d} t} \cdot \vec{A}(t,\vec{y}(t)) \right].
Taking the variation with respect to A^{\mu} you get Maxwell's equations with the current given by the single-particle current and variation with respect to \vec{y} gives you the equation of motion for a particle in this electromagnetic field.
The only trouble with this fully self-consistent equation is that it is plagued with the well-known problems of the self-consistent description of such a classical point-particle-field system. It's the so-called radiation-reaction problem. A very detailed analysis of the problem can be found in the book
F. Rohrlich, Classical Charged Particles, World Scientific (2007)
A good review is also given in Jackson's Classical Electrodynamics and in vol. 2 of Landau/Lifshitz.
Of course you can read off the usual approximations that can be solved in the sense of perturbation theory, i.e., you consider the motion of the particle in a given electromagnetic field or the electromagnetic field of a particle moving along a given trajectory \vec{y}(t), leading to the retarded potentials (Lienard-Wiechert potentials).
Another approach is to use kinetic theory or hydrodynamics (look, e.g., for "magneto hydrodynamics" in the usual textbooks).
