# Cellular Automata/Excitable media

## An introduction to excitable media

Excitable media are nonlinear dynamic systems known for exhibiting complex behavior that can be observed as pattern formation. They are usually defined by a reaction-diffusion differential equation.

${\displaystyle u_{t}(\mathbf {r} ,t)=D\nabla ^{2}u(\mathbf {r} ,t)+f(u(\mathbf {r} ,t))}$

The diffusion part provides stability and propagation of information, the reactive part provides interesting local dynamics.

A common example of excitable media are prey-predator systems. Such systems are described by a system of differential equations, one function for each of the observed protagonists.

${\displaystyle u_{t}(\mathbf {r} ,t)=D_{u}\nabla ^{2}u(\mathbf {r} ,t)+f(u,v)}$
${\displaystyle v_{t}(\mathbf {r} ,t)=D_{v}\nabla ^{2}v(\mathbf {r} ,t)+g(u,v)}$

We will discuss two different approaches to modeling excitable media. Discretization of differential equations and modeling with cellular automata.

### Boundary conditions

There are different ways to define boundary conditions for the reaction-diffusion equation.

Dirichlet boundary conditions

The value of the function at the boundary is given explicitly ${\displaystyle u(\mathbf {r_{0}} ,t)}$.

Cyclic boundaries

If the initial condition ${\displaystyle u(x,0)}$ is supposed to be periodic in space, cyclic boundary conditions can be used.

Zero-flux boundary conditions

If zero-flux is expected at the boundary than the component of the functions first derivative normal to the boundary is zero ${\displaystyle {\vec {n}}\cdot \nabla {u}=0}$ at the boundary. This can be achieved by reflecting function values from the inside over the boundary to the outside.

## Discretization of differential equations using the explicit FTCS method

The traditional method to simulate excitable media is discretization and numerical computation of the governing PDE. First the FTCS (forward-time centered-space method) discretization method is presented. Explicit methods are the simplest and the equations are similar to a cellular automaton, but are inadequate because of stability and convergence problems.

### Single PDE

We will first observe a single PDE describing a single function.

${\displaystyle u_{t}(\mathbf {r} ,t)=D\nabla ^{2}u(\mathbf {r} ,t)+f(u(\mathbf {r} ,t))}$

#### One-dimensional problem

In the one dimensional case the space vector becomes a single variable ${\displaystyle \mathbf {r} =x}$. The nabla operator becomes ${\displaystyle \nabla ^{2}={\frac {\partial ^{2}}{\partial {x^{2}}}}}$.

${\displaystyle {\frac {\partial u(x,t)}{\partial t}}=D{\frac {\partial ^{2}u(x,t)}{\partial x^{2}}}+f(u(x,t))}$

The partial differential equation is discretized.

${\displaystyle {\frac {u(x,t+\Delta {t})-u(x,t)}{\Delta {t}}}=}$
${\displaystyle \quad D{\frac {u(x+\Delta {x},t)-2u(x,t)+u(x-\Delta {x},t)}{\Delta {x}^{2}}}}$
${\displaystyle \quad +f(u(x,t))}$
Forward-time centered-space method

Each finite element at time ${\displaystyle t+\Delta {t}}$ is calculated from three neighboring elements at time ${\displaystyle t}$ (see figure at the right).

${\displaystyle u(x,t+\Delta {t})=u(x,t)+\,}$
${\displaystyle \quad +d(u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t))}$
${\displaystyle \quad +\Delta {t}f(u(x,t))}$

where the diffusion number is

${\displaystyle d=D{\frac {\Delta {t}}{\Delta {x}^{2}}}}$
Stability

The FTCS method is stable if

${\displaystyle d\leq 1/2}$
Boundary conditions

For periodic boundaries present values at the left boundary ${\displaystyle x=0}$ can be used to compute the future values at the right boundary ${\displaystyle x=L_{x}}$ and the other way round.

${\displaystyle u(0,t+\Delta {t})=u(0,t)+d(u(0+\Delta x,t)-2u(0,t)+u(L_{x},t))+\Delta {t}f(u(0,t))\,}$
${\displaystyle u(L_{x},t+\Delta {t})=u(L_{x},t)+d(u(0,t)-2u(L_{x},t)+u(L_{x}-\Delta x,t))+\Delta {t}f(u(L_{x},t))\,}$

If there is zero-flux ${\displaystyle \nabla u(0,t)=0}$ at the boundaries than values outside the boundary are reflections of values inside ${\displaystyle u(0+\Delta x,t)=u(0-\Delta x,t)}$.

${\displaystyle u(0,t+\Delta {t})=u(0,t)+2d(u(0+\Delta x,t)-u(0,t))+\Delta {t}f(u(0,t))\,}$
${\displaystyle u(L_{x},t+\Delta {t})=u(L_{x},t)+2d(u(L_{x},t)+u(L_{x}-\Delta x,t))+\Delta {t}f(u(L_{x},t))\,}$

#### Two-dimensional problem

In the two-dimensional case, the space vector becomes a variable pair ${\displaystyle \mathbf {r} =[x,y]}$. The nabla operator becomes ${\displaystyle \nabla ^{2}={\frac {\partial ^{2}}{\partial {x^{2}}}}+{\frac {\partial ^{2}}{\partial {y^{2}}}}}$.

${\displaystyle {\frac {\partial u(x,y,t)}{\partial t}}=D\left({\frac {\partial ^{2}u(x,y,t)}{\partial x^{2}}}+{\frac {\partial ^{2}u(x,y,t)}{\partial y^{2}}}\right)+f(u(x,y,t))}$

The partial differential equation is discretized using the forward-time centered-space method.

${\displaystyle {\frac {u(x,y,t+\Delta {t})-u(x,y,t)}{\Delta {t}}}=}$
${\displaystyle \quad =D{\frac {u(x+\Delta {x},y,t)-2u(x,y,t)+u(x-\Delta {x},y,t)}{\Delta {x}^{2}}}+}$
${\displaystyle \quad +D{\frac {u(x,y+\Delta {y},t)-2u(x,y,t)+u(x,y-\Delta {y},t)}{\Delta {y}^{2}}}+}$
${\displaystyle \quad +f(u(x,y,t))}$
Forward-time centered-space method

Each finite element at time ${\displaystyle t+\Delta {t}}$ is calculated from five neighboring elements at time ${\displaystyle t}$ (see figure at the right).

${\displaystyle u(x,y,t+\Delta {t})=u(x,t)+\,}$
${\displaystyle \quad +d_{x}(u(x+\Delta {x},y,t)-2u(x,y,t)+u(x-\Delta {x},y,t))+}$
${\displaystyle \quad +d_{y}(u(x,y+\Delta {y},t)-2u(x,y,t)+u(x,y-\Delta {y},t))+}$
${\displaystyle \quad +\Delta {t}f(u(x,y,t))}$

where the diffusion numbers are

${\displaystyle d_{x}=D{\frac {\Delta {t}}{\Delta {x}^{2}}}\quad d_{y}=D{\frac {\Delta {t}}{\Delta {y}^{2}}}}$
Stability

The FTCS method is stable if

${\displaystyle d_{x}\leq 1/2}$ and ${\displaystyle d_{y}\leq 1/2}$
Boundary conditions

The same ideas as in the one dimensional case can be used for two dimensions.

### System of PDE

A system of PDE describes two functions that interact with each other (prey-predator).

${\displaystyle u_{t}(\mathbf {r} ,t)=D_{u}\nabla ^{2}u(\mathbf {r} ,t)+f(u,v)}$
${\displaystyle v_{t}(\mathbf {r} ,t)=D_{v}\nabla ^{2}v(\mathbf {r} ,t)+g(u,v)}$

The interaction is local, which means, the dispersion part can be computed separately for each equation, and than the reaction part is added to the result.

## References

1. Joe D. Hoffman, Numerical Methods for Engineers and Scientists
2. Toffoli T., Margolus N., Cellular Automata Machines: A New Environment for Modeling, The MIT Press (1987), Cambridge, Massachusetts
3. Toffoli T., Cellular automata as an alternative to Differential equations, in Modeling Physics, Physica 10D, (1984)
4. http://www.jweimar.de/paper-abstracts.html
5. Robert Fisch, Janko Gravner, David Griffeath, Threshold-Range Scaling of Excitable Cellular Automata
6. Robert Fisch, Janko Gravner, David Griffeath, Metastability in the Greenberg-Hastings Model
7. Marcus R. Garvie Finite difference schemes for reaction-diffusion equations modeling predator-prey interactions in MATLAB