Cellular Automata/Examples of Plankton and Fish Dynamics

Definition of prey-predator systems
Prey-predator systems are a well known biological example of excitable media. This example is based upon an article (see references) about plankton and fish dynamics. The dynamics are given as a PDE system. The purpose of this example is to observe pattern formation and other dynamic properties of prey-predator systems using both finite difference equations and cellular automata.

PDE describing prey-predator systems
We define the prey-predator system using the a system of general reaction-diffusion differential equations.


 * $$ u_t(\mathbf{r},t) = D_u \nabla^2 u(\mathbf{r},t) + f(u,v) $$
 * $$ v_t(\mathbf{r},t) = D_v \nabla^2 v(\mathbf{r},t) + g(u,v) $$

The first equation describes the prey population $$u(\mathbf{r},t)$$ and the seccond equation describes the predator population $$v(\mathbf{r},t)$$. Nonlinear functions $$f(u,v)$$ and $$g(u,v)$$ describe local prey-predator dynamics, without them the system is decomposed into two separate diffusion equations.

We are interested in prey-predator dynamics of aquatic systems, where phytoplankton (small plants, food producers) is the prey and zooplankton (small animals) is the predator. Since we observe only slowly moving plankton the diffusion coefficient $$D_u=D_v=D$$ depends mainly on marine turbulence. The two nonlinear functions describing prey-predator dynamics are


 * $$ f(u,v) = P(u) - E(u,v) \, $$
 * $$ g(u,v) = \kappa E(u,v) -\mu v \, $$

$$P(u)$$ describes local growth and mortality of the prey, $$E(u,v)$$ describes the predation, $$\kappa$$ is food utilization and $$\mu$$ is the mortality rate of the predator. The selection of functions $$P(u)$$, $$E(u,v)$$ and parameters $$\kappa$$, $$\mu$$ depends on the observed population. Here a simple example is used, for a detailed description of the equations see the references.


 * $$ u_t(\mathbf{r},t) = D \nabla^2 u(\mathbf{r},t) + \frac{\alpha}{b}u(b-u) - \gamma\frac{u}{u+H}v $$
 * $$ v_t(\mathbf{r},t) = D \nabla^2 v(\mathbf{r},t) + \kappa\gamma\frac{u}{u+H}v - \mu v $$

$$\alpha$$ is the maximal growth rate of the prey, $$b$$ is the carrying capacity of the prey population, and H is the half-saturation abundance of prey.

Algebraic analysis
The first step is to introduce dimensionless variables


 * $$ \tilde{u}=\frac{u}{b} \quad \tilde{v}=\frac{v\gamma}{b\alpha} \quad \tilde{t}=\alpha{t} \quad \tilde{\mathbf{r}}=\mathbf{r} \left( \frac{\alpha}{D} \right)^{\frac{1}{2}} $$

and the new dimensionless parameters are


 * $$h=\frac{H}{b} \quad m=\frac{\mu}{\alpha} \quad k=\frac{\kappa\gamma}{\alpha}$$

the PDE system with only dimensionless quantities is (tildes are omitted)


 * $$ u_t = \nabla^2 u + u(1-u) - \frac{u}{u+h}v $$
 * $$ v_t = \nabla^2 v + k\frac{u}{u+h}v - mv $$

Local Dynamics


The local dynamics of the system are defined by the nonlinear reaction part of equations


 * $$ u_t = u(1-u) - \frac{u}{u+h}v \qquad  v_t = v + k\frac{u}{u+h}v - mv $$

We first observe stationary conditions $$ u_t = v_t = 0 \, $$ There are three stationary states $$(u,v)$$
 * total extinction $$(0,0)$$ is a sedal point for $$k,m,h>0$$
 * predator extinction $$(1,0)$$ is a sedal point if $$(u_*,v_*)$$ is in the biologically meaningful region $$u,v\geq 0$$ and a stable point otherwise
 * and coexistence $$(u_*,v_*)$$ can be a sedal or stationary point
 * $$ u_* = \frac{rh}{1-r} \quad r=m/k \quad v_* = (1-u_*)(h+u_*)  $$

Although the reaction dynamics depend on three parameters $$h$$ is the most important. If


 * $$ h < \frac{1-r}{r} $$

than the population is in a biologically meaningful region $$u_*,v_*\geq 0$$. If


 * $$ h < \frac{1-r}{1+r} $$

than there is a stable limit cycle surrounding the steady state. There is a figure of such a cycle where the path begins near from the stationary point $$(u_*,v_*+\Delta{v})$$ and the parameters are $$k=0.2$$, $$r=0.3$$ and $$h=0.4$$.

Numerical simulation using finite difference equations
A direct approach to simulate prey-predator dynamics in the plankton example is to discretize the reaction-diffusion differential equation, and to solve it numerically. The dynamics is first observed in an one dimensional space and than on two dimensions. In both cases using a zero-flux boundary $$\nabla u=\nabla v=0$$.

The PDE discretization method used is FTCS, which is the simplest. It is capable of producing complex behavior we would like to observe, but it has stability and convergence problems discussed below.

One dimension
The simulation is first performed in an one dimensional space of $$4000$$ elements. The discretization parameters are


 * $$ \Delta x = 1 \quad \Delta t = 0.2 \quad d=0.2 $$

The prey-predator system parameters are


 * $$ D=1 \quad k=2.0 \quad r=0.4 \quad h=0.3 $$

For initial conditions the uniform stationary point $$(u_*,v_*)$$ is used


 * $$ u(x,0) = u_* \, $$
 * $$ v(x,0) = v_* + \epsilon x + \delta \, $$

a gradual increase of predator population from the left to the right is added to ignite the spatial dynamics, the parameters are


 * $$ \epsilon=10^{-6} \quad \delta=-1.5 \cdot 10^{-3} $$

There are two figures showing the system at times 0, 640, 2600. At a certain point in time and space the species distribution becomes chaotic. The chaotic part shows a fast growth.



The system is not converging by decreasing the time difference $$\Delta{t}\rightarrow0$$. The source of the problem can be seen by observing the local dynamics of $$u(t)$$ and $$v(t)$$ inside the chaotic region $$x=1500$$ from time $$t=0$$ to $$t=640$$ using two different time differences $$\Delta{t}=0.20$$ and $$\Delta{t}=0.02$$. With decreasing the time difference the influence of diffusion grows over the influence of local dynamics. On the above figure the effect would be seen as a slower spreading of the chaotic region.



Two dimensions
The simulation is then performed in a two-dimensional space of $$576 \times 384$$ elements. The discretization parameters are


 * $$ \Delta x = \Delta y = 2 \quad \Delta t = 0.5 \quad d=0.125 $$

The prey-predator system parameters are


 * $$ D=1 \quad k=2.0 \quad r=0.4 \quad h=0.3 $$

For initial conditions the uniform stationary point $$(u_*,v_*)$$ is used


 * $$ u(x,0) = u_* + \epsilon_1(x-0.1y-144).*(x-0.1y-432) \, $$
 * $$ v(x,0) = v_* + \epsilon_2(x-450) - \epsilon_3(y-150) \, $$

a gradual increase of predator population is added to ignite the spatial dynamics, the parameters are


 * $$ \epsilon_1=2\cdot10^{-7} \quad \epsilon_2=3\cdot10^{-5} \quad \epsilon_3=1.2\cdot10^{-4} $$

2D dynamics are observed on time points $$t=\{0,100,150,200,300,400,1000\}$$. Spiral patterns can be seen arising from the initial gradient, at the end patterns brake and the system ends in chaos.



Another problem arises because the FTCS method is not unconditionally stable. Some results may become negative, and the functions describing local dynamics diverge very fast in the negative quadrants.