$\newcommand{\R}{\mathbb R}$You do not need to "discretize $\nabla$". Also, you wrote the diffusion equation incorrectly. The correct version is this: 
\begin{equation}
\frac{\partial f(r,t)}{\partial t}=\nabla\cdot[B(r,t)\,\nabla f(r,t)],
\end{equation}
where $f:=\phi$, $B:=D$, and $\cdot$ denotes the dot product. In the coordinate form, this equations is 
\begin{equation}
\frac{\partial f(r,t)}{\partial t}=\sum_{j=1}^n [B(r,t)\,(D_j^2 f)(r,t)+(D_j B)(r,t)\,
(D_j f)(r,t)],
\end{equation}
where $D_j$ is the operator of the partial differentiation with respect to the $j$th coordinate of $r\in\R^n$. 

Now discretization becomes straightforward, by replacing the partial derivatives by the corresponding differences: 
\begin{equation}
\frac{df_i(r,t)}{dt}=
\sum_j A_{i,j}[B_i(t)\,(f_j(t)-f_i(t))+
(B_j(t)-B_i(t))(f_j(t)-f_i(t))] 
\end{equation}
or, simply, 
\begin{equation}
\frac{df_i(r,t)}{dt}=
\sum_j A_{i,j}B_j(t)(f_j(t)-f_i(t)).  
\end{equation}

More generally, we can write 
\begin{equation}
\frac{df_i(r,t)}{dt}=
\sum_j A_{i,j}B_{i,j}(t)(f_j(t)-f_i(t)),  
\end{equation}
where the $B_{i,j}$'s are nonnegative functions. This will describe a general continuous-time random walk on the network. One may want to recall at this point that the diffusion equation describes an approximation of jump processes (which are continuous-time random walks) by processes continuous in time.