Barkley model
In the simplest case the equations for the Barkley model are
- <math>
- <math>
where the two reaction terms take the form
- <math>
- <math>
A more general form of the model is possible (see below). The role of the model parameters <math>\epsilon, a, b </math> is discussed in 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.
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>
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. The <math>v</math>-nullcline is also a straight line, but this is not essential to the model.
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.
<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.
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.
For the moment ignore diffusion and consider the ordinary differential equations for just the reaction terms in the model:
- <math>
- <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>
- <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>
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}
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>
Because <math>v</math> evolves on the slow time scale it may be time stepped by the explicit Euler method.
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>
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>
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>
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).
In the most general case the equations for the Barkley model are
- <math>
- <math>
- <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.)
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 here.
- Barkley D. (1991) "A model for fast computer simulation of waves in excitable media". Physica 49D, 61–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.
- 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.
Category:Cardiac Dynamics Category:Dynamical Systems Category:Eponymous __AUTOLINKER{1|evolution|dynamics|language}