Barkley model
| Dwight Barkley (2008), Scholarpedia, 3(11):1877. | doi:10.4249/scholarpedia.1877 | revision #91030 [link to/cite this article] |
The Barkley model 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.
In the simplest case the equations for the Barkley model are
\[ \frac{\partial u }{ \partial t} = f(u,v) + \nabla^2 u, \] \[ \frac{\partial v }{ \partial t} = g(u,v), \]
where the two reaction terms take the form
\[ f(u,v) = \frac{1}{\epsilon}u (1-u)\left(u-\frac{v+b}{a}\right), \] \[ g(u,v) = u - v. \]
A more general form of the model is possible (see below). The role of the model parameters \(\epsilon, a, b \) is discussed in section 1.1.
Figure 1 and Figure 2 show typical waves in two and three dimensions from simulations of the model.
Contents |
Nullclines
The nullclines for the Barkley model are displayed in Figure 3 where one sees the defining feature of the model: the nullclines for the nonlinear \(u\) reaction kinetics are straight lines. The \(u\)-nullclines are given by \(f(u,v)=0\) so that the three branches are
\[ u = \begin{cases} 0, \\ (v + b)/a, \\ 1 \end{cases} \]
The middle branch (dashed in Figure 3) sets the excitation threshold. In practice only the N-shaped portion shown in bold in Figure 3 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 \(f(u,v)\) is a cubic polynomial in \(u\) and hence is everywhere smooth. The nullclines in Figure 3 should be compared with the S-shaped nullclines for the FitzHugh-Nagumo model. The \(v\)-nullcline is also a straight line, but this is not essential to the model.
Parameters
The parameter \(\epsilon\) sets the timescale separation between the fast \(u\)-equation and the slow \(v\)-equation, and is thus taken to be small. Larger \(a\) gives longer excitation duration and increasing \(b/a\) gives a larger excitability threshold.
The spiral waves in Figure 1 are for parameter values \(a = 0.75\ ,\) \(b=0.02\ ,\) \(\epsilon = 0.02\ .\) The simulation domain is of size \(L = 150\) on a side. Colors indicate the values of the \(v\)-field. The scroll waves in Figure 2 are for parameter values \(a = 0.8\ ,\) \(b=0.01\ ,\) \(\epsilon = 0.02\ .\) The simulation domain is of size \(L = 30\) on a side. Shown are isosurfaces corresponding to \(u = 1/2\ .\) Figure 4 shows the dynamics of a single spiral wave as a function of parameters \(a\) and \(b\) for fixed \(\epsilon=0.02\ .\) The spirals may undergo period rotations or various types of meander.
Why are simulations fast?
Figure 5 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 \(u\) sampled along a one-dimensional cut through a rotating spiral wave. Plotted on the left are values in the \((u,v)\) plane along this cut. Points indicate the numerical grid used in the simulation. Blue signifies that \(u < 10^{-3}\ .\) 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 \(u\)-nullcline. All other dynamics are slow by comparison. The simplicity of the \(u\)-nullclines permits efficiently, accurately, and stably within these fast regions.
- Large quiescent regions. There are substantial quiescent regions, shown in blue, in which \(u\simeq 0\ .\) In these regions the right-hand-side of the \(u\) 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 Figure 6. 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: \[ \frac{d u }{d t} = \frac{1}{\epsilon} u (1-u)\left(u-\frac{v+b}{a}\right), \] \[ \frac{d v }{d t} = u - v. \]
Letting \(u^n\) and \(v^n\) denote the values of \(u\) and \(v\) at time step \(n\ ,\) one explicit Euler time step of the reaction kinetics is given by
\[ u^{n+1} = u^n + \frac{\triangle t}{\epsilon} \, u^n (1-u^n)(u^n - u_{th}), \] \[ v^{n+1} = v^n + \triangle t \, (u^n - v^n), \]
where \(u_{th} = (v^n + b)/a\) and \(\triangle t\) 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:
\[ 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} \]
where \(u\) 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 \(u^{n+1}\) gives
\[\tag{1} 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} \]
By expanding the denominators in the above expressions it can be seen that
this scheme agrees with the explicit-Euler for small \(\triangle t\)
and has the same order of accuracy. However, unlike the explicit form, the
semi-implicit form is numerically stable for arbitrarily large
\(\triangle t\) and in the limit of large \(\triangle t\)
it goes over to:
\[ u^{n+1} = \begin{cases} 0 & \mbox {if } u^n < u_{th} \\ 1 & \mbox {if } u^n > u_{th} \end{cases} \] i.e. essentially a 2 state cellular-automaton.
Because \(v\) evolves on the slow time scale it may be time stepped by the explicit Euler method.
Efficient treatment of diffusion
The \(u\)-field is approximately constant (and zero) over substantial regions of the computational domain, e.g. Figure 5. With the approximation that \(u=0\) for \(u<\delta\ ,\) where \(\delta\) is a small numerical parameter, the Laplacian of the \(u\)-field is zero in these regions and need not be evaluated. To avoid unnecessary computation, the Laplacian of \(u\) 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:
\[ 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} \]
where \(i,j\) labels grid points and \(h\) is the grid spacing.
A scatter evaluation of the sums \(\Sigma_{ij}\) at all grid points is obtained by looping over the grid indices and scattering values of \(u\) to neighboring points:
\[ \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} \]
where all \(\Sigma_{ij}\) are initially zero. \(\Sigma_{ij}/h^2\) gives an approximation to the Laplacian at point \((i,j)\ .\)
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, \(\Sigma^0_{ij}\) and \(\Sigma^1_{ij}\ ,\) 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:
\[ \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} \]
where \(F(u_{ij}, u_{th})\) is given by equation (1)
and \(\nu \equiv \triangle t/h^2\ ,\) and where \(r\) has a
value of 0 or 1, and \(s\) 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
\[ \frac{\partial u }{ \partial t} = f(u,v) + \nabla^2 u, \] \[ \frac{\partial v }{ \partial t} = g(u,v) + D_v \nabla^2 v, \] where \[ f(u,v) = \frac{h(v)}{\epsilon}u (1-u)(u-u_{th}(v)). \]
\(D_v\) 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 \(g(u,v)\ ,\) \(h(v)\ ,\) and \(u_{th}(v)\) can be specified as needed. The original and simplest choice is \(h(v)=1\ ,\) \(u_{th}(v)=(v+b)/a\ ,\) and \(g(u,v)=u-v\ .\) Another important case suggested by Bär and Eiswirth (1993) is the original choice of \(h\) and \(u_{th}\) but with \(g(u,v)\) an nonlinear function of \(u\ .\) For example with \(g(u,v) = u^3 - v\) the model can exhibit spiral breakup as shown in Figure 7. (Parameters \(a = 0.75\ ,\) \(b=0.0006\ ,\) and \(\epsilon\) 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 here.
References
- 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.
Internal references
- 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.
