Cell invasion (4) : Numerical solution of drift-diffusion equation

Back

The cell invasion problem requires an accurate numerical solution of coupled partial differential equations. The fields (concentration profiles) f(z,t) will depend on one spatial coordinate z and the time t. The spatial coordinates will be restricted to the intervall [0,len], with well-defined boundary conditions at the left and right limits: The current must dissapear outside the intervall.

I will try to solve this problem by discretizing the z-coordinate into a regular grid,

z_k = k\;\Delta z,

and then treating each z_k as a separate dynamical variable. The problem can thereby be mapped onto a set of coupled ordinary differential equations, which are solved with the Runga-Kutta method.

As a simple test and practice case, consider the drift-diffusion equation

\frac{\partial}{\partial t}c(z,t) = - v \frac{\partial}{\partial z}c(z,t) + D \frac{\partial^2}{\partial z^2}c(z,t),

with position- and time-independent diffusion constant D and drift velocity v.

For the following tests, we set the initial distribution c(z,t) equal to a narrow, normalized Gaussian (variance=10), placed at the center of the intervall. The spatial resolution is set to dz=1.

Test 1: Free diffusion in a box

First we set the drift-velocity equal zero. The initial peak should broaden with time. Asymptotically, a flat distribution should develop, while the integral (total particle number) is preserved. The simulation shows the expected behaviour:



Test 2: Diffusion with drift (boundaries far away)

The case of combined diffusion and drift has been compared directly with the analytical solution (solid lines). One expects the center z_c(t) of the Gaussian to move according to

z_c(t) = v\;t

and the variance to grow linearly

\sigma^2(t)=\sigma^2(0)+ 2Dt

The simulation (dashed lines) fits the analytical curves nicely, even for the modest spatial resolution of dz=1:



Test 3: Diffusion with drift against a wall

Finally, we choose a negative drift velocity and wait until the distribution is pushed againts the left wall. We now expect the stationary solution

c(z)=\frac{v}{D}\;e^{-(v/D)z}.

Again, the simulation (symbols) agrees very well with the analytical solution (solid lines):



As a conclusion, the simple numerical method seems to work as expected. I have also checked that the results do not change significantly when the resolution dz is reduced.

Next