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

Local dynamics cycle (starting in (u*,v*))

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.

Spatial distributions of populations in the 1D case at time t=640

Spatial distributions of populations in the 1D case at time t=2640

The system is not converging by decreasing the time difference $\Delta {t}\rightarrow 0$. 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.

Local dynamics (u(t),v(t)) observed at x=1500 t=0-640 with Δt=0.20

Local dynamics (u(t),v(t)) observed at x=1500 t=0-640 with Δt=0.02

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\cdot 10^{-7}\quad \epsilon _{2}=3\cdot 10^{-5}\quad \epsilon _{3}=1.2\cdot 10^{-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.

t=0

t=100

t=150

t=200

t=300

t=400

t=1000

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.

References

1. Alexander B. Medvinsky, Sergei V. Petrovskii, Irene A. Tikhonova, Horst Malchow, Bai-Lian Li, Spatiotemporal Complexity of Plankton and Fish Dynamics, SIAM Review, Volume 44, Number 3, pp. 311-370

1. Marcus R. Garvie, Finite difference schemes for reaction-diffusion equations modeling predator-prey interactions in MATLAB PDF