An introduction to excitable media
edit
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 .
u
t
(
r
,
t
)
=
D
∇
2
u
(
r
,
t
)
+
f
(
u
(
r
,
t
)
)
{\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.
u
t
(
r
,
t
)
=
D
u
∇
2
u
(
r
,
t
)
+
f
(
u
,
v
)
{\displaystyle u_{t}(\mathbf {r} ,t)=D_{u}\nabla ^{2}u(\mathbf {r} ,t)+f(u,v)}
v
t
(
r
,
t
)
=
D
v
∇
2
v
(
r
,
t
)
+
g
(
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
edit
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
u
(
r
0
,
t
)
{\displaystyle u(\mathbf {r_{0}} ,t)}
.
Cyclic boundaries
If the initial condition
u
(
x
,
0
)
{\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
n
→
⋅
∇
u
=
0
{\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
edit
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
edit
We will first observe a single PDE describing a single function.
u
t
(
r
,
t
)
=
D
∇
2
u
(
r
,
t
)
+
f
(
u
(
r
,
t
)
)
{\displaystyle u_{t}(\mathbf {r} ,t)=D\nabla ^{2}u(\mathbf {r} ,t)+f(u(\mathbf {r} ,t))}
One-dimensional problem
edit
In the one dimensional case the space vector becomes a single variable
r
=
x
{\displaystyle \mathbf {r} =x}
. The nabla operator becomes
∇
2
=
∂
2
∂
x
2
{\displaystyle \nabla ^{2}={\frac {\partial ^{2}}{\partial {x^{2}}}}}
.
∂
u
(
x
,
t
)
∂
t
=
D
∂
2
u
(
x
,
t
)
∂
x
2
+
f
(
u
(
x
,
t
)
)
{\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.
u
(
x
,
t
+
Δ
t
)
−
u
(
x
,
t
)
Δ
t
=
{\displaystyle {\frac {u(x,t+\Delta {t})-u(x,t)}{\Delta {t}}}=}
D
u
(
x
+
Δ
x
,
t
)
−
2
u
(
x
,
t
)
+
u
(
x
−
Δ
x
,
t
)
Δ
x
2
{\displaystyle \quad D{\frac {u(x+\Delta {x},t)-2u(x,t)+u(x-\Delta {x},t)}{\Delta {x}^{2}}}}
+
f
(
u
(
x
,
t
)
)
{\displaystyle \quad +f(u(x,t))}
Forward-time centered-space method
Each finite element at time
t
+
Δ
t
{\displaystyle t+\Delta {t}}
is calculated from three neighboring elements at time
t
{\displaystyle t}
(see figure at the right).
u
(
x
,
t
+
Δ
t
)
=
u
(
x
,
t
)
+
{\displaystyle u(x,t+\Delta {t})=u(x,t)+\,}
+
d
(
u
(
x
+
Δ
x
,
t
)
−
2
u
(
x
,
t
)
+
u
(
x
−
Δ
x
,
t
)
)
{\displaystyle \quad +d(u(x+\Delta x,t)-2u(x,t)+u(x-\Delta x,t))}
+
Δ
t
f
(
u
(
x
,
t
)
)
{\displaystyle \quad +\Delta {t}f(u(x,t))}
where the diffusion number is
d
=
D
Δ
t
Δ
x
2
{\displaystyle d=D{\frac {\Delta {t}}{\Delta {x}^{2}}}}
Stability
The FTCS method is stable if
d
≤
1
/
2
{\displaystyle d\leq 1/2}
Boundary conditions
For periodic boundaries present values at the left boundary
x
=
0
{\displaystyle x=0}
can be used to compute the future values at the right boundary
x
=
L
x
{\displaystyle x=L_{x}}
and the other way round.
u
(
0
,
t
+
Δ
t
)
=
u
(
0
,
t
)
+
d
(
u
(
0
+
Δ
x
,
t
)
−
2
u
(
0
,
t
)
+
u
(
L
x
,
t
)
)
+
Δ
t
f
(
u
(
0
,
t
)
)
{\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))\,}
u
(
L
x
,
t
+
Δ
t
)
=
u
(
L
x
,
t
)
+
d
(
u
(
0
,
t
)
−
2
u
(
L
x
,
t
)
+
u
(
L
x
−
Δ
x
,
t
)
)
+
Δ
t
f
(
u
(
L
x
,
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
∇
u
(
0
,
t
)
=
0
{\displaystyle \nabla u(0,t)=0}
at the boundaries than values outside the boundary are reflections of values inside
u
(
0
+
Δ
x
,
t
)
=
u
(
0
−
Δ
x
,
t
)
{\displaystyle u(0+\Delta x,t)=u(0-\Delta x,t)}
.
u
(
0
,
t
+
Δ
t
)
=
u
(
0
,
t
)
+
2
d
(
u
(
0
+
Δ
x
,
t
)
−
u
(
0
,
t
)
)
+
Δ
t
f
(
u
(
0
,
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))\,}
u
(
L
x
,
t
+
Δ
t
)
=
u
(
L
x
,
t
)
+
2
d
(
u
(
L
x
,
t
)
+
u
(
L
x
−
Δ
x
,
t
)
)
+
Δ
t
f
(
u
(
L
x
,
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
edit
In the two-dimensional case, the space vector becomes a variable pair
r
=
[
x
,
y
]
{\displaystyle \mathbf {r} =[x,y]}
. The nabla operator becomes
∇
2
=
∂
2
∂
x
2
+
∂
2
∂
y
2
{\displaystyle \nabla ^{2}={\frac {\partial ^{2}}{\partial {x^{2}}}}+{\frac {\partial ^{2}}{\partial {y^{2}}}}}
.
∂
u
(
x
,
y
,
t
)
∂
t
=
D
(
∂
2
u
(
x
,
y
,
t
)
∂
x
2
+
∂
2
u
(
x
,
y
,
t
)
∂
y
2
)
+
f
(
u
(
x
,
y
,
t
)
)
{\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 .
u
(
x
,
y
,
t
+
Δ
t
)
−
u
(
x
,
y
,
t
)
Δ
t
=
{\displaystyle {\frac {u(x,y,t+\Delta {t})-u(x,y,t)}{\Delta {t}}}=}
=
D
u
(
x
+
Δ
x
,
y
,
t
)
−
2
u
(
x
,
y
,
t
)
+
u
(
x
−
Δ
x
,
y
,
t
)
Δ
x
2
+
{\displaystyle \quad =D{\frac {u(x+\Delta {x},y,t)-2u(x,y,t)+u(x-\Delta {x},y,t)}{\Delta {x}^{2}}}+}
+
D
u
(
x
,
y
+
Δ
y
,
t
)
−
2
u
(
x
,
y
,
t
)
+
u
(
x
,
y
−
Δ
y
,
t
)
Δ
y
2
+
{\displaystyle \quad +D{\frac {u(x,y+\Delta {y},t)-2u(x,y,t)+u(x,y-\Delta {y},t)}{\Delta {y}^{2}}}+}
+
f
(
u
(
x
,
y
,
t
)
)
{\displaystyle \quad +f(u(x,y,t))}
Forward-time centered-space method
Each finite element at time
t
+
Δ
t
{\displaystyle t+\Delta {t}}
is calculated from five neighboring elements at time
t
{\displaystyle t}
(see figure at the right).
u
(
x
,
y
,
t
+
Δ
t
)
=
u
(
x
,
t
)
+
{\displaystyle u(x,y,t+\Delta {t})=u(x,t)+\,}
+
d
x
(
u
(
x
+
Δ
x
,
y
,
t
)
−
2
u
(
x
,
y
,
t
)
+
u
(
x
−
Δ
x
,
y
,
t
)
)
+
{\displaystyle \quad +d_{x}(u(x+\Delta {x},y,t)-2u(x,y,t)+u(x-\Delta {x},y,t))+}
+
d
y
(
u
(
x
,
y
+
Δ
y
,
t
)
−
2
u
(
x
,
y
,
t
)
+
u
(
x
,
y
−
Δ
y
,
t
)
)
+
{\displaystyle \quad +d_{y}(u(x,y+\Delta {y},t)-2u(x,y,t)+u(x,y-\Delta {y},t))+}
+
Δ
t
f
(
u
(
x
,
y
,
t
)
)
{\displaystyle \quad +\Delta {t}f(u(x,y,t))}
where the diffusion numbers are
d
x
=
D
Δ
t
Δ
x
2
d
y
=
D
Δ
t
Δ
y
2
{\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
d
x
≤
1
/
2
{\displaystyle d_{x}\leq 1/2}
and
d
y
≤
1
/
2
{\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
edit
A system of PDE describes two functions that interact with each other (prey-predator).
u
t
(
r
,
t
)
=
D
u
∇
2
u
(
r
,
t
)
+
f
(
u
,
v
)
{\displaystyle u_{t}(\mathbf {r} ,t)=D_{u}\nabla ^{2}u(\mathbf {r} ,t)+f(u,v)}
v
t
(
r
,
t
)
=
D
v
∇
2
v
(
r
,
t
)
+
g
(
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.
Other PDE discretization methods
edit
Modeling with cellular automata
edit
Greenberg-Hastings Model
edit
Conclusions
edit
References
edit
Joe D. Hoffman, Numerical Methods for Engineers and Scientists
Toffoli T., Margolus N., Cellular Automata Machines: A New Environment for Modeling , The MIT Press (1987), Cambridge, Massachusetts
Toffoli T., Cellular automata as an alternative to Differential equations, in Modeling Physics, Physica 10D, (1984)
http://www.jweimar.de/paper-abstracts.html
Robert Fisch, Janko Gravner, David Griffeath, Threshold-Range Scaling of Excitable Cellular Automata
Robert Fisch, Janko Gravner, David Griffeath, Metastability in the Greenberg-Hastings Model
Marcus R. Garvie Finite difference schemes for reaction-diffusion equations modeling predator-prey interactions in MATLAB