Given a point $\vec{r}$, divide the space in two pieces: 1) $V$ such that $\color{#0000ff}{\Large\vec{r} \in V}$ and 2) $\tilde{V}$. $V \cup \tilde{V}$ is the whole space. At any time $t$, the number of particles $N_{V}\left(t\right)$ in $V$ plus the number of particles $N_{\tilde{V}}\left(t\right)$ in $\tilde{V}$ is a constant of motion: $N_{V}\left(t\right) + N_{\tilde{V}}\left(t\right) = \mbox{constant}\,, \forall\ t$, such that
$$
{{\rm d}N_{V}\left(t\right) \over {\rm d}t}
+
{{\rm d}N_{\tilde{V}}\left(t\right) \over {\rm d}t} = 0\,,
\qquad
\forall\ t
$$
$N_{V}\left(t\right)$ is given by
$$
N_{V}\left(t\right) = \int_{V}\rho\left(\vec{R},t\right)\,{\rm d}V
\qquad\Longrightarrow\qquad
{{\rm d}N_{V}\left(t\right) \over {\rm d}t}
=
{\LARGE\int}_{\!\!\!\!\!\!V}\!\!\!\!\!\!{\partial\rho\left(\vec{R},t\right) \over \partial t}\,{\rm d}V
$$
where $\rho\left(\vec{R},t\right)$ is, at time $t$, the particle density at point $\vec{R}$. In an infinitesimal ${\rm d}S$ surface element, $\mbox{at}\ \vec{R}$, of $V$ and during the time interval $\Delta t$, the number of particles that cross ${\rm d}S$ is given by
$\Delta N_{\tilde{V}}\left(t\right)
\equiv
\rho\left(\vec{R}, t\right)\left\{\left[\vec{v}\left(\vec{R}, t\right)
\Delta t\right]{\rm d}S\cos\left(\angle\right)\right\}\quad$. $\vec{v}\left(\vec{R}, t\right)$ is the flux velocity and $\angle$ is the angle between $\vec{v}\left(\vec{R}, t\right)$ and a normalized vector $\vec{n}$ perpendicular to ${\rm d}S$. $\left\vert\vec{n}\right\vert = 1$. Then
${\rm d}N_{\tilde{V}}\left(t\right)
=
\vec{\rm J}\left(\vec{R}, t\right)\cdot{\rm d}\vec{S}\Delta t$
where
$\vec{\rm J}\left(\vec{R}, t\right)
\equiv
\rho\left(\vec{R}, t\right)\vec{v}\left(\vec{R}, t\right)$ and
${\rm d}\vec{S} \equiv {\rm d}S\,\,\vec{n}\quad$.
$\vec{\rm J}\left(\vec{R}, t\right)$ is the $Current\ Density\ Vector$.
$$
{{\rm d}N_{\tilde{V}}\left(t\right) \over {\rm d}t}
=
\lim_{\Delta t \to 0}\int_{S}{\Delta N_{\tilde{V}}\left(t\right) \over \Delta t}
=
\int_{S}\vec{\rm J}\left(\vec{R}, t\right)\cdot{\rm d}\vec{S}
=
\int_{V}\nabla\cdot\vec{\rm J}\left(\vec{R}, t\right)\,{\rm d}V
$$
where we have used the Gauss Divergence Theorem. Then
$$
{\LARGE\int}_{\!\!\!\!\!\!V}\!\left[%
{\partial\rho\left(\vec{R},t\right) \over \partial t}
+
\nabla\cdot\vec{\rm J}\left(\vec{R}, t\right)
\right]\,{\rm d}V
=
0
$$
\begin{align}
&\lim_{V \to 0}{1 \over V}\!\!\!{\LARGE\int}_{\!\!\!\!\!\!V}\!\!\!\!\left[%
{\partial\rho\left(\vec{R},t\right) \over \partial t}
+
\nabla\cdot\vec{\rm J}\left(\vec{R}, t\right)
\right]\,{\rm d}V
=
0\
\Longrightarrow\
\color{#ff0000}{\large%
{\partial\rho\left(\vec{r},t\right) \over \partial t}
+
\nabla\cdot\vec{\rm J}\left(\vec{r}, t\right)
=
0\,\color{#0000ff}{;\forall\, \vec{r}, t}}
\end{align}