View source for Barkley model

[[Image:spiral_waves.png|thumb|200px|right|F1|Spiral waves from the
Barkley model.]]
The <strong>Barkley model</strong> is a system of reaction-diffusion equations modeling
excitable media and oscillatory media. The model is often used as a qualitative model in pattern forming systems like the Belousov-Zhabotinsky reaction and other systems that are well described by the interaction of an activator and an inhibitor component. It is very similar to the
FitzHugh-Nagumo (FHN) model and for many purposes the Barkley model can be
considered to be the same as the FHN model.  The Barkley model and associated
numerical algorithm were originally proposed as a method for performing very
fast numerical simulations of spiral waves in two dimensions and scroll waves
in three dimensions. Numerical simulations of the model remain fast to this
day, but, in addition, the model has gained importance because many studies of
waves in excitable media have been based on the Barkley model and its
extensions.

[[Image:scroll_wave.png|thumb|200px|right|F2|Scroll waves from the
Barkley model.]] 

In the simplest case the equations for the Barkley model are

:<math>
\frac{\partial u }{ \partial t} = f(u,v) +  \nabla^2 u, 
</math>
:<math>
\frac{\partial v }{ \partial t} = g(u,v),
</math>

where the two reaction terms take the form

:<math>
f(u,v) = \frac{1}{\epsilon}u (1-u)\left(u-\frac{v+b}{a}\right),
</math>
:<math>
g(u,v) = u - v.
</math>

A more general form of the model is possible (see 
[[#General form of the model|below]]).  
The role of the model parameters <math>\epsilon, a, b </math> is discussed 
in [[#Parameters|section 1.1]].

<figref>Spiral_waves.png</figref> and <figref>Scroll_wave.png</figref> show typical waves in two and three
dimensions from simulations of the model.

== Nullclines ==

[[Image:nullclines.png|thumb|200px|left|F3|Nullclines for the Barkley
model. Bold blue shows the ''N''-shaped <math>u</math>-nullcline and green
shows the <math>v</math>-nullcline. ]]

The nullclines for the Barkley model are displayed in  <figref>Nullclines.png</figref>
where one sees the defining feature of the model: the nullclines for the
nonlinear <math>u</math> reaction kinetics are straight lines.  The
<math>u</math>-nullclines are given by <math>f(u,v)=0</math> so that the three
branches are

:<math>
u = 
\begin{cases} 
  0,  \\
  (v + b)/a, \\
  1
\end{cases}
</math>

The middle branch (dashed in  <figref>Nullclines.png</figref>) sets the excitation
threshold.  In practice only the ''N''-shaped portion shown in bold in 
<figref>Nullclines.png</figref> is relevant, since for spiral and scroll wave solutions the
system does not pass through the corners where branches of the nullclines
intersect.  It should be emphasized that while the nullclines are piecewise
linear, the function <math>f(u,v)</math> is a cubic polynomial in
<math>u</math> and hence is everywhere smooth.  The nullclines in 
<figref>Nullclines.png</figref> should be compared with the ''S''-shaped nullclines for the
[[FitzHugh-Nagumo model|FitzHugh-Nagumo model]].
The <math>v</math>-nullcline is also a straight line, but this is not
essential to the model.

[[Image:Barkley_model_flowers.png|right|thumb|200px|F4|Top: yellow region denotes periodically rotating spirals and cyan denotes meandering spirals. Bottom: Cut through the parameter space at <math>b=0.05</math> with the different states illustrated by tip paths.]]
=== Parameters ===

The parameter <math>\epsilon</math> sets the timescale separation between the
fast <math>u</math>-equation and the slow <math>v</math>-equation, and is thus
taken to be small.  Larger <math>a</math> gives longer excitation duration and
increasing <math>b/a</math> gives a larger excitability threshold.

The spiral waves in  <figref>Spiral_waves.png</figref> are for parameter values
<math>a = 0.75\ ,</math> <math>b=0.02\ ,</math> <math>\epsilon = 0.02\ .</math> The
simulation domain is of size <math>L = 150</math> on a side. Colors indicate the
values of the <math>v</math>-field. 
The scroll waves in  <figref>Scroll_wave.png</figref> are for parameter values
<math>a = 0.8\ ,</math> <math>b=0.01\ ,</math> <math>\epsilon = 0.02\ .</math> The
simulation domain is of size <math>L = 30</math> on a side. Shown are 
isosurfaces corresponding to <math>u = 1/2\ .</math>
<figref>Barkley_model_flowers.png</figref> shows the dynamics of a single spiral wave as a function of parameters <math>a</math> and <math>b</math> for fixed <math>\epsilon=0.02\ .</math> The spirals may undergo period rotations or various types of 
[[Meander of spiral waves|meander]].

== Why are simulations fast? ==

<figref>Barkley_model_cut_fine.gif</figref> illustrates the typical behavior of variables in
the Barkley model as a spiral or scroll wave evolves.  Plotted on the right is
the fast variable <math>u</math> sampled along a one-dimensional cut through a
rotating spiral wave.  Plotted on the left are values in the
<math>(u,v)</math> plane along this cut.  Points indicate the numerical grid
used in the simulation.  Blue signifies that <math>u < 10^{-3}\ .</math>  Red
highlights values at a fixed grid point.

[[Image:barkley_model_cut_fine.gif|right|thumb|400px|F5|1D cut through
well-resolved simulation of spiral waves. The grid spacing is <math>\triangle
x=0.25</math> and the time step is <math>\triangle t=0.018\ .</math> Each frame
of the animation corresponds to 5 time steps.]]
[[Image:barkley_model_cut_coarse.gif|right|thumb|400px|F6|1D cut through a
simulation of spiral waves with large space and time steps.  The
grid spacing is <math>\triangle x=0.75</math> and the time step is
<math>\triangle t=0.16\ .</math> Each frame of the animation corresponds to one
time step.  This simulation is nearly 100 times faster than that in 
<figref>Barkley_model_cut_fine.gif</figref>.]]

From these plots one can see the two factors exploited in simulations of the
Barkley model:

*''Time-scale separation''.  The variable u at a given grid point evolves on a fast time scale only when a wave front or wave back passes a grid point. This fast evolution corresponds to moving between the two vertical branches of the <math>u</math>-nullcline.  All other dynamics are slow by comparison. The simplicity of the <math>u</math>-nullclines permits efficiently, accurately, and stably within these fast regions.

*''Large quiescent regions''.  There are substantial quiescent regions, shown in blue, in which <math>u\simeq 0\ .</math> In these regions the right-hand-side of the <math>u</math> equation effectively vanishes and neither diffusion nor the nonlinear reaction kinetics need to be evaluated, thereby reducing the number of computations which need to be performed each time step. 

Combining these feature in a single algorithm means that floating-point
computations can be minimized and simulations can be run stably at large
space and time steps as in  <figref>Barkley_model_cut_coarse.gif</figref>.  For large time steps the
algorithm effectively reduces to a cellular-automaton type simulation of
excitable media in which the fast variable 
takes on just the values 0 and 1 (Gerhardt, Schuster, and Tyson; 1990).

== Algorithm ==

=== Semi-implicit time stepping of the reaction kinetics ===

For the moment ignore diffusion and consider the ordinary differential
equations for just the reaction terms in the model:
:<math>
\frac{d u }{d t} = \frac{1}{\epsilon} u (1-u)\left(u-\frac{v+b}{a}\right),
</math>
:<math>
\frac{d v }{d t} = u - v.
</math>

Letting <math>u^n</math> and <math>v^n</math> denote the values of 
<math>u</math> and <math>v</math> at time step <math>n\ ,</math> 
one explicit Euler time step of the reaction kinetics is given by

:<math>
u^{n+1} = u^n + \frac{\triangle t}{\epsilon} \, u^n (1-u^n)(u^n - u_{th}),
</math>
:<math>
v^{n+1} = v^n + \triangle t \, (u^n - v^n),
</math>

where <math>u_{th} = (v^n + b)/a</math> and <math>\triangle t</math> is the
time step. 

With the explicit-Euler method, the time step is limited by a numerical
stability constraint.  In order to take larger time steps, a semi-implicit
approximation may be used. The simplest form is:

:<math>
u^{n+1} = 
\begin{cases}
u^n + (\triangle t/\epsilon) \, u^{n+1} (1-u^n)(u^n - u_{th})  
& \mbox{if } u^n \le u_{th} \\
u^n + (\triangle t/\epsilon) \, u^n (1-u^{n+1})(u^n - u_{th})  
& \mbox{if } u^n > u_{th} 
\end{cases}
</math>

where <math>u</math> at the future time is used in those factors on the
right-hand side that undergo largest relative change over the time
step. Solving for <math>u^{n+1}</math> gives

:<math>\label{Fdef}
u^{n+1} = F(u^n, u_{th}) = 
\begin{cases}
u^n/\left (1 - (\triangle t/\epsilon) \, (1-u^n)(u^n - u_{th}) \right ) 
& \mbox {if } u^n \le u_{th} \\
\left (u^n + (\triangle t/\epsilon) \, u^n (u^n - u_{th}) \right )/
\left (1 + (\triangle t/\epsilon) \, u^n (u^n - u_{th}) \right )
& \mbox {if } u^n > u_{th} 
\end{cases}
</math>


By expanding the denominators in the above expressions it can be seen that
this scheme agrees with the explicit-Euler for small <math>\triangle t</math>
and has the same order of accuracy.  However, unlike the explicit form, the
semi-implicit form is numerically stable for arbitrarily large 
<math>\triangle t</math> and in the limit of large <math>\triangle t</math> 
it goes over to:

:<math>
u^{n+1} = 
\begin{cases}
0 & \mbox {if } u^n < u_{th} \\
1 & \mbox {if } u^n > u_{th} 
\end{cases}
</math>
i.e. essentially a 2 state cellular-automaton. 

Because <math>v</math> evolves on the slow time scale it may be time stepped
by the explicit Euler method.

=== Efficient treatment of diffusion ===

The <math>u</math>-field is approximately constant (and zero) over substantial
regions of the computational domain, e.g.  <figref>Barkley_model_cut_fine.gif</figref>.  With the
approximation that <math>u=0</math> for <math>u<\delta\ ,</math> where
<math>\delta</math> is a small numerical parameter, the Laplacian of the
<math>u</math>-field is zero in these regions and need not be evaluated. To
avoid unnecessary computation, the Laplacian of <math>u</math> can be
evaluated by scattering values rather than by gathering values as follows.

Consider the five-point finite-difference formula for the Laplacian on a
regular square lattice in 2D:

:<math>
h^2 \nabla^2 u_{ij} \simeq \Sigma_{ij} 
\equiv u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 4 u_{i,j} 
</math>

where <math>i,j</math> labels grid points and <math>h</math> is the grid
spacing. 

A scatter evaluation of the sums <math>\Sigma_{ij}</math> at all grid points
is obtained by looping over the grid indices and scattering values of
<math>u</math> to neighboring points:

:<math>
\begin{array}{l}
\mbox{For Each}\ i, j \\
~~~~~ \Sigma_{ij} = \Sigma_{ij} - 4 u_{ij}  \\
~~~~~ \Sigma_{i-1,j} = \Sigma_{i-1j} + u_{ij}  \\
~~~~~ \Sigma_{i+1,j} = \Sigma_{i+1j} + u_{ij}  \\
~~~~~ \Sigma_{i,j-1} = \Sigma_{ij-1} + u_{ij}  \\
~~~~~ \Sigma_{i,j+1} = \Sigma_{ij+1} + u_{ij}  
\end{array}
</math>

where all <math>\Sigma_{ij}</math> are initially zero.  
<math>\Sigma_{ij}/h^2</math> gives an approximation to the
Laplacian at point <math>(i,j)\ .</math>

What makes a scatter evaluation of the sum desirable is that it can be
combined into an algorithm for the reaction dynamics in such a way that
unnecessary computation is avoided at points which make no contribution to the
Laplacian of the u-field.  Specifically, two arrays,
<math>\Sigma^0_{ij}</math> and <math>\Sigma^1_{ij}\ ,</math> are employed.  One
is used to perform the current update while the other is used to collect the
sum of neighboring points for use at the next time step.  The following
algorithm advances the solution on a two-dimensional grid forward one time
step, while simultaneously preparing for the next time step:

:<math>
\begin{array}{l}
\mbox{For Each}\ i, j \\
~~~~~\mbox{If}\ u_{ij} < \delta~ \mbox{Then} \\
~~~~~ ~~~~~ u_{ij} = \nu \, \Sigma^r_{ij} \\
~~~~~ ~~~~~ v_{ij} = (1 - \triangle t) v_{ij} \\
~~~~~\mbox{Else} \\
~~~~~ ~~~~~ u_{th} = (v_{ij} + b)/a \\
~~~~~ ~~~~~ v_{ij} = v_{ij} + \triangle t (u_{ij} - v_{ij}) \\
~~~~~ ~~~~~ u_{ij} = F(u_{ij}, u_{th}) + \nu \, \Sigma^r_{ij} \\
~~~~~ ~~~~~ \Sigma^{s}_{ij} = \Sigma^{s}_{ij} - 4 u_{ij}  \\
~~~~~ ~~~~~ \Sigma^{s}_{i-1,j} = \Sigma^{s}_{i-1,j} + u_{ij} \\
~~~~~ ~~~~~ \Sigma^{s}_{i+1,j} = \Sigma^{s}_{i+1,j} + u_{ij} \\
~~~~~ ~~~~~ \Sigma^{s}_{i,j-1} = \Sigma^{s}_{i,j-1} + u_{ij} \\
~~~~~ ~~~~~ \Sigma^{s}_{i,j+1} = \Sigma^{s}_{i,j+1} + u_{ij} \\
~~~~~\Sigma^{r}_{ij} = 0 \\
\mbox{Swap values of}\ r \mbox{ and}\ s
\end{array}
</math>


where <math>F(u_{ij}, u_{th})</math> is given by equation \eqref{Fdef}
and <math>\nu \equiv \triangle t/h^2\ ,</math> and where <math>r</math> has a
value of 0 or 1, and <math>s</math> has the other value. These values are
swapped at the end of each time step.

This simple algorithm illustrates the essence of the Barkley model without
taking into account boundary or initial conditions. More advanced algorithms
are available and are recommended (Dowle, Mantel, and Barkley; 1997).

== General form of the model ==

In the most general case the equations for the Barkley model are

[[Image:Spiral_breakup.png|thumb|right|F7|Spiral-wave breakup in the
Barkley model.]]

:<math>
\frac{\partial u }{ \partial t} = f(u,v) +  \nabla^2 u, 
</math>
:<math>
\frac{\partial v }{ \partial t} = g(u,v) + D_v \nabla^2 v, 
</math>
where 
:<math>
f(u,v) = \frac{h(v)}{\epsilon}u (1-u)(u-u_{th}(v)).
</math>

<math>D_v</math> is the diffusion constant for the slow species, or
equivalently the ratio of diffusion coefficients since the model is written in
space units in which the diffusion coefficient of the fast species is one.
The three functions <math>g(u,v)\ ,</math> <math>h(v)\ ,</math> and
<math>u_{th}(v)</math> can be specified as needed. The original and simplest
choice is <math>h(v)=1\ ,</math> <math>u_{th}(v)=(v+b)/a\ ,</math> and
<math>g(u,v)=u-v\ .</math> Another important case suggested by Bär and Eiswirth 
(1993) is the original choice of
<math>h</math> and <math>u_{th}</math> but with <math>g(u,v)</math> an
nonlinear function of <math>u\ .</math>  For example with <math>g(u,v) = u^3 -
v</math> the model can exhibit spiral breakup as shown in 
<figref>Spiral_breakup.png</figref>. (Parameters <math>a = 0.75\ ,</math> <math>b=0.0006\ ,</math> and
<math>\epsilon</math> changed from 0.06 to 0.08.)

== Software ==

Two open source computer programs EZ-Spiral and EZ-Scroll implement the
Barkley model in two and three dimensions, respectively. The programs are written in the C programming language and
employ OpenGL graphics. All spiral and scroll wave images appearing in this
article were generated using this software. These programs can be found
[http://www.warwick.ac.uk/staff/D.Barkley/Software/ez_software.html here].

== References == 

* Barkley D. (1991) "A model for fast computer simulation of waves in excitable media". ''Physica'' '''49D''', 61&ndash;70.
* Dowle M., Mantel R.M., and Barkley D. (1997) "Fast simulations of waves in three-dimensional excitable media". ''Int. J. Bif. Chaos'' '''7''', 2529--2545.
* Gerhardt M., Schuster H., and Tyson J.J. (1990) "A Cellular Automaton Model of Excitable Media Including Curvature and Dispersion". ''Science'' '''247''', 1563.   
* M. Bär and M. Eiswirth (1993) "Turbulence due to spiral breakup in a continuous excitable medium". ''Phys. Rev. E'' '''48''', R1635 -- R1637.

<!-- Authors, please check this list and remove any references that are irrelevant. This list is generated automatically to reflect the links from your article to other accepted articles in Scholarpedia. -->
<b>Internal references</b>

* Anatol M. Zhabotinsky (2007) [[Belousov-Zhabotinsky reaction]]. Scholarpedia, 2(9):1435.

* Richard J. Field (2008) [[Chemical reaction kinetics]]. Scholarpedia, 3(10):4051.

* Vladimir S. Zykov (2008) [[Excitable media]]. Scholarpedia, 3(5):1834.

* Eugene M. Izhikevich and Richard FitzHugh (2006) [[FitzHugh-Nagumo model]]. Scholarpedia, 1(9):1349.

* Gregoire Nicolis and Anne De Wit (2007) [[Reaction-diffusion systems]]. Scholarpedia, 2(9):1475.

* Philip Holmes and Eric T. Shea-Brown (2006) [[Stability]]. Scholarpedia, 1(10):1838.


==See also==

[[Category:Cardiac Dynamics]]
[[Category:Dynamical Systems]]
[[Category:Eponymous]]
__AUTOLINKER{1|evolution|dynamics|language}