Cellular Automata/Examples of Plankton and Fish Dynamics

From Wikibooks, open books for an open world
< Cellular Automata
Jump to: navigation, search

Definition of prey-predator systems[edit]

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[edit]

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[edit]

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[edit]

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[edit]

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[edit]

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}\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.

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[edit]

The simulation is than performed in an 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.

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[edit]

  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