Well, the von Neumann equation holds within Schroedinger's picture and it is immediate to prove in quantum physics (differently from the Liouville equation which needs to preventively establish the non-trivial Liouville theorem for the symplectic measure on the space of phases).
Indeed, as $\rho$ is an incoherent superposition of pure states, one has,
$$\rho = \sum_k p_k |\psi_k \rangle \langle \psi_k|\:.$$
Hence, $\rho$ evolves in time as effect of the standard Schroedinger evolution of pure states of the mixture,
$$\rho(t) = \sum_k p_k U_t |\psi_k \rangle \langle \psi_k| U^\dagger_t = U_t \rho U^\dagger_t \:,$$
where $U_t = e^{-\frac{it}{\hbar}H}$ is the usual time evolutor. This identity immediately leads to von Neumann equation,
$$\frac{d}{dt} \rho(t) = -\frac{i}{\hbar}[H, \rho]\:,\tag{1}$$
where the derivative is computed with respect to the strong operator topology (or the weak one), the same notion of time derivative as that used in Schroedinger equation. Unfortunately, and in my opinion erroneously, that derivative is very often indicated by $\frac{\partial}{\partial t}$ instead of $\frac{d}{dt}$. This would make sense if $\rho$ included another time dependence. Think of Heisenberg picture, when operators already depend on time in Schroedinger picture. In that case to distinguish between $\partial/\partial t$ (referred to some parametric time dependence already present in Schroedinger picture) and $d/dt$ (referred to the total time dipendence, including the one arising from Heisenberg evolution) makes sense. Here, instead, we have only one time dependence.
If $\rho(t)$ were a time-depending observable in Schroedinger picture, passing in Heisenberg picture we would obtain
$$\frac{d\rho_H}{dt} = \frac{\partial \rho_H}{\partial t} +\frac{i}{\hbar}[H, \rho_H]=0\:,\tag{2}$$
i.e.,
$$\frac{d\rho_H}{dt}=0\:.$$
Actually all that is trivial without computing any derivative, since
$$\rho_H(t) = U^\dagger_t \rho(t) U_t = U^\dagger_t U_t \rho U^\dagger_t U_t= \rho\:.$$
The conserved quantity would be $\rho$ itself, but it does not make much physical sense a priori, because $\rho$ is a state, not an observable and Heisenberg evolution is not appropriate for it. (Even if the result is formally correct since states, in Heisenberg picture, are constant in time.)
In classical physics the partial derivative $\partial/\partial t$ always makes sense, since $\rho$ is a function of both $t$ and the state in the space of phases. In that formulation (2) has a precise meaning in terms of a conservation law. It says that, along the integral lines of Hamiltonian flow, the density of probability is constant. In other words, the Hamiltonian flow is incompressible (see the pair of final remarks in my answer to this question Conserved quantities and total derivatives?)
OP's question (v1) is essentially asking
Does the operator identity
$$ e^{\frac{it}{\hbar}[\hat{H},~\cdot~]}\hat{A}~ =~ e^{i\hat{H}t/\hbar}\hat{A}e^{-i\hat{H}t/\hbar} \tag{1} $$
have an analog using functions/symbols $H$ and $A$ rather than operators $\hat{H}$ and $\hat{A}$, respectively?
The answer is: Yes, in terms of the Groenewold-Moyal star product $\star$. If the Poisson bracket is written as
$$ \{A,B\}_{PB}~:=~A \stackrel{\leftarrow}{\partial_I} \omega^{IJ} \stackrel{\rightarrow}{\partial_J} B, \qquad \partial_I~\equiv~\frac{\partial}{\partial z^I}, \qquad \{z^I,z^J\}_{PB}~=~\omega^{IJ}, \tag{2} $$
where $z^1, \ldots, z^{2n}$ are canonical coordinates, and $A,B$ are functions/symbols (as opposed to operators), then the star product reads
$$ A\star B~:=~A \exp\left(\stackrel{\leftarrow}{\partial_I}\frac{i\hbar}{2} \omega^{IJ} \stackrel{\rightarrow}{\partial_J}\right) B~=~AB+\frac{i\hbar}{2}\{A,B\}_{PB} +{\cal O}(\hbar^2). \tag{3} $$
And then an analog of eq. (1) is
$$ e^{\frac{it}{\hbar}[~H~\stackrel{\star}{,}~\cdot~]}A~=~ e_{\star}^{iHt/\hbar}\star A\star e_{\star}^{-iHt/\hbar}, \tag{4} $$
where
$$[A\stackrel{\star}{,}B]~:= A\star B-B\star A~=~i\hbar\{A,B\}_{PB} +{\cal O}(\hbar^3)\tag{5} $$
is the star commutator, and
$$e_{\star}^B~:=~1+B\sum_{n=1}^{\infty}(\star B)^{n-1}/n!\tag{6} $$
is the star exponential.
Best Answer
If the Hamiltonian is time-dependent the evolution of a pure state is $$ |\psi(t)\rangle = {\mathcal T} exp\left[ -\frac{i}{\hbar} \int_{t_0}^t d\tau H(\tau)\right] |\psi(0)\rangle = U(t; t_0) |\psi(0)\rangle $$ where $\mathcal T$ is the time-ordering operator and $$ i\hbar\dot{U} = H(t) U(t, t_0) $$ Then a density matrix evolves according to $$ \rho(t) = U(t; t_0) \rho(0) U^\dagger (t; t_0) $$ and you can check that taking the time derivative gives $$ i\hbar\dot{\rho} = \left[ H(t), \rho(t)\right] $$