The large core limit of spiral waves in excitable media: A numerical approach

The large core limit of spiral waves in excitable media: A numerical approach

Sebastian Hermann School of Mathematics and Statistics, University of Sydney, NSW 2006, Australia (shermann@maths.usyd.edu.au), supported by an Endeavour Australia-Europe Award, an EIPRS award and the Deutscher Akademischer Austauschdienst.    Georg A. Gottwald School of Mathematics and Statistics, University of Sydney, NSW 2006, Australia (georg.gottwald@sydney.edu.au), supported by the Australian Research Council.
Abstract

We modify the freezing method introduced by Beyn & Thümmler, 2004, for analyzing rigidly rotating spiral waves in excitable media. The proposed method is designed to stably determine the rotation frequency and the core radius of rotating spirals, as well as the approximate shape of spiral waves in unbounded domains. In particular, we introduce spiral wave boundary conditions based on geometric approximations of spiral wave solutions by Archimedean spirals and by involutes of circles. We further propose a simple implementation of boundary conditions for the case when the inhibitor is non-diffusive, a case which had previously caused spurious oscillations.

We then utilize the method to numerically analyze the large core limit. The proposed method allows us to investigate the case close to criticality where spiral waves acquire infinite core radius and zero rotation frequency , before they begin to develop into retracting fingers. We confirm the linear scaling regime of a drift bifurcation for the rotation frequency and the core radius of spiral wave solutions close to criticality. This regime is unattainable with conventional numerical methods.

e

xcitable media, pattern formation, equivariance, unbounded domains

{AMS}

35B36, 65M99, 35K57

1 Introduction

Excitable media are abundant in nature, and appear in physical, chemical and biological systems. Prominent examples include cAMP waves in slime mold aggregation [39] and intracellular calcium waves [4], electrical waves in cardiac and nerve tissue [46, 9], and the auto-catalytic Belousov-Zhabotinsky reaction [45]. Excitable media involve the interwoven dynamics of so called activators and inhibitors. An important class of excitable media contains non-diffusive inhibitors. This is the classic situation in cardiac dynamics where the inhibitor consists of (relatively) immobile ion channels. Excitable media support localized pulses and periodic wave trains. In two dimensions rotating vortices (or spirals) are possible [45], and in three dimensions scroll waves occur [45, 47, 49, 30, 31].

In this work we focus on rigidly rotating spiral waves of single-diffusive excitable media. These spiral waves perform a time-periodic motion where the spiral tip moves on a perfect circle around an unexcited region, the spiral core. Such spirals are characterized by their rotation frequency and their core radius . Varying the excitability of the medium leads to qualitative changes in the type of motion the spiral is performing. For certain parameters rigidly rotating spirals bifurcate via a Hopf bifurcation into meandering spirals [48, 3, 22]. A different type of bifurcation occurs at low excitabilities where rigidly rotating spirals develop into travelling fingers with an overall zero curvature. Approaching this bifurcation, the rotation frequency of a rigidly rotating spiral wave decreases down to zero and its core radius diverges. Past this bifurcation rigidly rotating spirals do not exist. Whereas the bifurcation to meandering spirals is well-understood numerically and theoretically [3, 51, 14, 17, 37, 38], the latter one is not. Theoretically, several attempts have been made to study this large core limit, using kinematic theory [53, 43, 12, 32, 22, 23, 52, 13], dynamical systems theory [1] and non-perturbative asymptotics [20]. Computational results, on the other hand, are rare. In [20] the large core limit was investigated in terms of the growing velocity of the spiral wave tip close to criticality. Using a non-perturbative approach it was shown that the growing velocity behaves linearly as a function of the excitability close to criticality. This was verified numerically.

Here we set out to tackle the numerically much more difficult problem of determining the rotation frequency and the radius of the core of a spiral wave in this limit. This presents a huge computational challenge. Given current computer power, determining the rotation frequency and the core radius from direct simulations of the underlying reaction-diffusion system is restricted to small core spirals. In the large core limit, where the rotation frequency becomes zero and the core radius becomes infinite, no reliable results have been obtained so far.

Our aim here is twofold. First we present a numerical method to study rigidly rotating spirals in excitable media which is able to capture the spiral wave even in the large core limit without having to solve for computationally expensive large domains. We then apply this method to obtain the rotation frequency and the core radius of spiral wave solutions in the large core limit, and establish their scaling behaviour close to criticality to study the type of bifurcation.

The numerical method we present is based on the freezing method [6] in which the dynamics of an equivariant partial differential equation is split orthogonally into the shape dynamics and the dynamics of the associated symmetry group. This method has been successfully applied to many types of equivariant systems [5, 40, 42]. However, when applying the freezing method to single-diffusive excitable media several problems were encountered [6, 41, 8]. Firstly, the inhibitor exhibits spurious oscillations, the amplitude of which increases when approaching the large core limit. It was suggested that these oscillations are caused by numerical instabilities linked to the mixed hyperbolic-parabolic nature of the problem. The oscillations could be partially controlled by an upwind-downwind scheme for small core spirals. Secondly, the application of Neumann boundary conditions causes the shape of the spiral wave to deviate near the boundary from the form expected in an unbounded domain. However, it turns out that although the shape is not properly resolved near the boundary, the group parameters, i.e. rotation frequency and core radius, are accurately reproduced for small core spirals. This can be understood on a phenomenological basis by considering that the spiral wave coils shield the solution near the core, and information only flows outwards for rotating spiral wave solutions.
We will present an implementation of the freezing method which will overcome these difficulties. At first, we will identify two separate types of spatial oscillations, which are caused by different mechanisms: oscillations in the interior of the computational domain linked to numerical instabilities of the non-diffusive inhibitor, and oscillations close to the boundary caused by the application of Neumann boundary conditions. We eliminate the former by employing a semi-implicit Crank-Nicolson scheme and the latter by imposing a subtle implementation of a different type of boundary condition. We implement so-called spiral wave boundary conditions, approximating the spiral wave by an Archimedean spiral or by an involute of a circle, which approximately respect the symmetry of the solution, and produce the correct shape of a spiral wave in an effective unbounded domain at the boundary. We find that the spatial oscillations can be controlled if (a) boundary conditions which respect the symmetry of the solution are employed, and (b) the implementation of these boundary conditions is such that the boundary is coupled to the interior.

The paper is organized as follows. In Section 2 we present the excitable media model under consideration. In Section 3 we describe the original freezing method. This method is then modified to suit excitable media with a non-diffusive inhibitor in Section 4, where we introduce spiral wave boundary conditions and their discrete implementations. In Section 5 we apply our method to study the large core limit and determine the scaling behaviour of the rotation frequency and the core radius of spiral wave solutions close to criticality for varying excitability . We conclude with a discussion in Section 6.

2 Model

We consider here the Barkley model [2] for an activator and an inhibitor described by

(1)
(2)

Although the numerical method we will describe in Sections 3 and 4 is independent of the particular model used, we illustrate some basic properties of excitable media with the Barkley model (1). Our choice of model is motivated by the fact that it incorporates the ingredients of an excitable system in a compact and lucid way. Thus, for the rest state is linearly stable with decay rates along the activator direction and along the inhibitor direction. Perturbing above the threshold (in 0D) will lead to growth of . In the absence of the inhibitor the activator will saturate at leading to a bistable system. The positive inhibitor growth factor forces the activator to decay back to . Finally also the inhibitor with the refractory time constant will decay back to . For with the system is in 0D no longer excitable but instead bistable.

We have included a diffusion term for the inhibitor in (1). However, we will be concentrating on the case of vanishing diffusivity for with . This case is more relevant for cardiac dynamics where the inhibitor models the (relatively) immobile potassium and sodium ion channels. Moreover, we will see in Section 5 that problems of applying the freezing method as proposed in [6] to excitable media are caused by the lack of coupling of the boundary and the interior when .

We shall fix in our numerical simulations , and and vary the excitability parameter if not stated otherwise.

The Barkley model supports, in a well-defined parameter region [2], rigidly rotating spirals. These spiral wave solutions are characterized by their rotation frequency and their core radius . In the large core limit the rotation frequency approaches zero and the core radius becomes infinite at a critical value of the excitability parameter . At criticality a rigidly rotating spiral wave becomes a finger propagating in the transverse direction only, with the speed of the corresponding travelling wave. Finger like initial conditions will curl up and eventually develop into rigidly rotating spirals below criticality, or into retracting fingers above criticality (see Figure 1).

(a) Finger developing into a spiral
(b) Retracting finger
Figure 1: Temporal evolution of a retracting and spiraling finger solution of the Barkley model (1) at 3 different time steps. The motion is clockwise for the spiraling finger and from right to left for the retracting finger. The white curve represents the trace of the tip as defined by (3).

Whereas the rotation frequency is well defined for rigidly rotating spirals, the definition of the core radius is less clear (see [21] for a recent discussion). We define the core of a rigidly rotating spiral to be the maximal region which has never been excited during the course of one revolution of the spiral. To estimate this region we describe the boundary of the core as the trace of the tip during one revolution where the tip is defined as the intersection of the two contour lines

(3)

of the activator and the inhibitor, respectively [2]. The value solves with . This definition ensures that the distance of the tip to the centre of the spiral is close to minimal. Note that this definition is arbitrary. In the large core limit, however, the value of the core radius calculated as the radius of the circle traced by the spiral tip does not vary much proportionally for different choices of the contour lines used to define the location of the tip.

3 Freezing method

In this Section we briefly outline the so called freezing method as introduced in [6] before we propose our modifications to this method in Section 4.

We consider reaction-diffusion systems on the plane,

(4)

with , , , and a diffusion matrix with constant coefficients. The system (4) is equivariant under the action of the special Euclidean group consisting of rotations and translations. The Barkley model (1) described in the previous section belongs to this class of equations with .

Equivariant systems of the form (4) can be solved by the freezing method introduced in [6] for equivariant systems. This beautiful numerical method utilizes the fact that the dynamics of equivariant systems can be decomposed into two parts: the dynamics on the symmetry group and the dynamics orthogonal to it. Equivariant systems can be cast into a skew product form whereby the dynamics on the symmetry group is driven by the so called shape dynamics which is orthogonal to the group dynamics. This idea was developed in the seminal paper [3] and then put on a more rigorous footing in [28, 14, 15, 17, 37, 36, 38, 18, 35]. The method can be applied for any symmetry group. Here we restrict the exposition of the method to equivariance with respect to the Euclidean group, and follow closely [6, 41]. For a more general framework of the method we refer to [5, 7, 40, 42].

3.1 General Setup

Let be an element of the special Euclidean group consisting of the angle of rotation and translation . Two group elements are linked by the operation with

The action of the Euclidean group on functions can be defined by

The group action satisfies the properties

where is the unit element of and the identity matrix. Note that system (4) is equivariant under the action of , i.e.

(5)

The invariance of equation (4) with respect to the Euclidean group implies that it is possible to construct new solutions from a given solution by applying symmetry operations. In particular, we can rewrite a solution as

(6)

with . By formally differentiating (6) with respect to we obtain, upon using the equivariance condition (5),

which can be rearranged to yield

(7)

where we used and set . The derivative of the group action with respect to can be calculated as

where and . The gradient acts on vector-valued functions as . This yields the expression

(8)

Introducing new group variables by setting and defining parameters and allows for an elimination of the group variable in (8). Using this transformation we write (8) as

(9)

The equations for the group variables are given by

(10)
(11)

For constant equation (11) describes a rotation on a circle with centre at

(12)

and radius of rotation

(13)

for some point on the circle. The core radius can be determined by applying equation (13) to the tip of the spiral (as defined by (3)), using the centre of the core given by (12). The core radius can thus be calculated from the group parameters, which are calculated in the freezing procedure.

So far, the path on the group in (7) is arbitrary. Therefore three additional degrees of freedom, equaling the dimension of the group , exist. To close the system we fix the location on the group orbit, and augment the equations by a phase condition

(14)

where the functional maps into . The phase condition can be chosen in several ways. For the time–dependent problem (7) we will use the condition

(15)

which assigns the location on the group orbit by minimizing the temporal change of at each time step. Note that this condition is equivalent to requiring that is orthogonal to the group orbit at . This condition allows us to determine the freezing parameters , which are implicitly contained in (15) through . These values are then subsequently fed into the shape dynamics to update (see Section 3.2.1 for details on the implementation).

When considering the stationary problem of (7), we will use the following slightly simpler phase condition

(16)

which determines the location on the group orbit by minimizing the distance of to some template function . The phase condition (16) is independent of , and can therefore not be used directly in the time-dependent setting with a semi-implicit Crank-Nicolson scheme. Note that in order to use the template function must be sufficiently close to a solution of (4). For further details on phase conditions see [7, 16].

We summarize the closed system of partial differential algebraic equations (7), (9) – (11) and (14). The group dynamics is given by

(17)
(18)

which is driven through the parameters by the shape dynamics

(19)
(20)

with given by (9). This system is called the frozen system since for rigidly rotating spirals its solution evolves into a stationary solution for .

Our aim is to numerically determine the shape of rigidly rotating spiral wave solutions of the Barkley model (1) and the corresponding core radius and rotation frequency. Note that the shape dynamics (19) entirely determines the solution and the parameters which can then be used to determine the rotation frequency and the centre of the core via (12). To this end, it would be sufficient to solve the stationary problem corresponding to (19), using for example a Newton solver (as done in [42]). However, the high-dimensionality of the stationary problem, needed for an accurate resolution of the solution, requires a sufficiently good guess for the initialization of the Newton-Raphson method which otherwise would not converge. The initial guess for the Newton solver will be generated by the application of the freezing method for the time-dependent problem (19).

3.2 Discretization

We implement the freezing method for the Barkley model (1) with for Cartesian coordinates and for polar coordinates. The choice of the coordinate system depends on the parameter range and the specific situation. Cartesian coordinates are better suited for the investigation of the large core limit where a finger-like solution can intersect the boundary almost perpendicularly and Neumann boundary conditions are a good approximation for the unbounded domain. Polar coordinates are better suited for small core radii where several revolutions of the spiral are usually inside the computational domain.

Simulations in Cartesian coordinates are performed on a rectangular domain with points and gridsizes and . The shape dynamics reads for Cartesian coordinates as

(21)
(22)
(23)

with . The group variables can be determined by equations (17).

Simulations in polar coordinates with are performed on a rectangular domain consisting of grid points with gridsizes and . The analogue to (21) reads as

(24)
(25)
(26)

and the dynamics of the group variables is again given by (17).

The temporal integration for both (21) and (24), is done in discrete time steps . We denote the values of the fields and at the -th time-step and spatial location (or ) with , by and , respectively. Spatial derivatives are evaluated using second-order central differences.

3.2.1 Time-dependent freezing

We use a second order semi-implicit Crank-Nicolson scheme to solve the time-dependent problem (21) (or (24)) whereby the linear terms are treated implicitly and the nonlinear term is treated explicitly with an Adams-Bashforth scheme [34, 11]. The nonlinear freezing term can be rendered as an effectively linear term to be included into the Crank-Nicolson part of the temporal discretization, by first obtaining through the phase condition

The additional computational costs of solving the semi-implicit equations is far outweighed by the less restrictive time-step required for the semi-implicit method when compared to an explicit Euler method.

In [6] an explicit Euler scheme with Neumann boundary conditions was used which exhibited spurious spatial oscillations in both, the interior and at the boundary of the domain. An upwind-downwind scheme was introduced to control the oscillations with partial success for large excitabilities at small values of only. The semi-implicit scheme does not exhibit oscillations in the interior of the domain. However, spatial oscillations at the boundary persist. In Sections 4.2 and 4.3 we will present an explanation for the oscillations at the boundary, and provide a simple method to eliminate them.

3.2.2 Stationary freezing

Travelling waves and spiral waves are both examples of relative equilibria [25, 19]. Relative equilibria have the representation with a time-independent . Hence the temporal evolution of a solution can be described entirely by the dynamics on the group. We may therefore find the spiral wave solutions of (21) by solving its associated stationary problem

(27)

or, analogously for polar coordinates and (24)

(28)

with . We solve the system for the unknowns . Typically we use and for Cartesian coordinates and and for polar coordinates. To solve this high-dimensional nonlinear system we use the Newton-Raphson method with line searches and backtracking (see for example [34]) with a terminating condition .

4 Boundary Conditions

In the introduction, two issues were described which, so far, prevented a successful stable application of the freezing method to excitable media. Firstly, numerical instabilities in the form of spatial oscillations spoiled results in the case of a non-diffusive inhibitor. In the previous Section we have eliminated spatial oscillations in the interior of the domain by using a semi-implicit Crank-Nicolson scheme in the time-dependent freezing problem (21) (or (24)). Spatial oscillations near the boundaries, however, persist. Similarly, solving the stationary freezing problem (27) (or (28)) leads to spatial oscillations near the boundary. Secondly, the usual Neumann boundary conditions do not reproduce the correct shape of a spiral in an unbounded domain. These problems are associated with the type of boundary condition used, and how they are formulated in the discretization.

Naturally, computations are performed on a bounded domain and appropriate boundary conditions have to be chosen. The boundary conditions should preferably reflect the nature of the investigated system and its solutions. We are interested here in approximating spiral wave solutions in unbounded domains. By unbounded domains we mean either the infinite limit, where spiral wave solutions do not decay at infinity, or finite domains with boundary conditions but where the computational domain is much smaller than the actual physical domain and the physical boundaries can be ignored. In these cases the usual Neumann boundary conditions are, in general, not well suited. We discuss their impact on the freezing method in Section 4.1. In Section 4.2 we present spiral boundary conditions and show how they can be implemented for freezing methods. The standard implementation of both these boundary conditions leads to oscillations of the inhibitor near the boundary. In Section 4.3 we will therefore introduce an implementation for both Neumann and spiral wave boundary conditions, which does not exhibit spurious spatial oscillations.

4.1 Neumann boundary conditions

Most work on spiral waves in excitable media uses either Dirichlet or Neumann boundary conditions (NBC) (e.g. [2]). These boundary conditions are physically meaningful for simulations of excitable media on bounded domains, e.g. in chemical experiments. However, if simulations are performed with the intention to understand the behaviour of spirals in unbounded domains, they have the disadvantage of not respecting the underlying symmetry at the boundary. Nevertheless, Neumann boundary conditions have been used extensively. For simplicity we restrict our discussion to polar coordinates, for which Neumann boundary conditions are formulated as and at . These are discretized according to

(29)

for . The values of can then be used in the evaluation of the diffusion and advection terms on the boundary.

(a) NBC
(b) SBC (35) for and
(c) SBC (35) for , and (49) for
Figure 2: Inhibitor of a frozen spiral wave solution calculated via the stationary problem (28). From left to right: representation in the Cartesian plane , in the polar -plane, and a close-up of the radial boundary. (a)–(c): Neumann boundary condition (NBC) (29), (d)–(f): spiral boundary condition (35) for both, activator and inhibitor, (g)–(i): spiral boundary condition (SBC) (35) for activator only and free boundary condition (49) for the inhibitor. We chose in the Barkley model (1), and , and for the spatial discretization.

In the two top left panels of Figure 2 we show contour plots of the inhibitor calculated by solving the stationary problem (28) in polar coordinates. Figure 2(a) shows the spiral solution on the circular domain with radius in the Cartesian plane . In Figure 2 we show the same solution in the -plane. One sees clearly how the contours bend towards the boundary to satisfy the Neumann boundary condition , at . This kink is localized near the boundary and does not extend far into the domain. However, the bending of the spiral wave solution leads to an effective wider transversal cross-section of the solution near the boundary. A wider activator profile allows the inhibitor to adopt larger amplitudes (cf. the darker color (online red) of the inhibitor in the close–up in Figure 2). The boundaries, however, do not affect the dynamics near the core – provided they are located sufficiently far away from the spiral wave tip – and hence the presence of a kink in the solution near the boundary does not affect the values of the rotation frequency and the core radius .

More important, from a numerical stability point of view, is the following issue. In the close–up of the inhibitor near the boundary in 2, strong spatial oscillations are clearly visible. For a more detailed view we show a slice of the inhibitor near the boundary along the grid line in Figure 3. These oscillations are a well known problem of the freezing method for excitable media, and are associated with the non-diffusive nature of the inhibitor [6, 41, 8]. These spatial oscillations occur only in the non-diffusive inhibitor , and are absent for the diffusive activator . We have checked that there are no oscillations for sufficiently large diffusion coefficient of the inhibitor.

The amplitude and the extent of the oscillations increase strongly for larger values of , which prohibits an accurate numerical analysis of the large core limit.

Figure 3: Left: Slice of the close–up in Figure 2 along a radial grid line for the inhibitor solution where NBCs are used for activator and inhibitor , demonstrating spurious spatial oscillations near the boundary. Right: The same slice but now with SBC for the activator and the one-sided derivative (49) as boundary condition for the inhibitor , corresponding to Figure 2.

4.2 Spiral wave boundary conditions

In this Section we introduce two boundary conditions based on simple geometrical approximations of spiral waves in order to mimic the behaviour of spirals in unbounded domains, or in finite domains when the size of the computational domain is much smaller than the physical domain. A natural version of this type of boundary conditions using Archimedean spirals was introduced in [10]. Therein, spiral wave solutions were “grown” with a predefined wavelength . We will expand and generalize this idea with the aim to apply it to the freezing method. Here the wavelength is not known a priori and, moreover, the centre of the spiral wave generally moves during the freezing procedure. In addition to the approximation by an Archimedean spiral, we will present boundary conditions where the spiral wave is approximated by an involute of a circle [45]. The two approximations coincide in the far field but differ near the spiral wave tip. We will formulate spiral boundary conditions for polar and Cartesian coordinates. It has to be noted that both, the Archimedean spiral and the involute of a circle, are just geometrical approximations of contour lines of the actual spiral wave solution in unbounded domains. There exists, to our knowledge, no rigorous theoretical justification for these approximations. However, we will see that these approximations may serve as convenient constructs to formulate boundary conditions.

We assume a spiral wave solution with a constant far field wavelength , centred at the origin of the polar coordinate system . A common choice to describe such a spiral wave geometrically is by means of an Archimedean spiral. In this approximation contour lines of the spiral wave are given by

(30)

allowing us to simplify . The parameter is given by , which can easily be seen by noting that a spiral wave profile is -periodic with , and therefore . The approximation of a spiral wave solution of an excitable medium by an Archimedean spiral is fairly accurate in the far field, away from the spiral wave core, but fails close to it. This is illustrated in Figure 4, where two examples of rigidly rotating spiral wave solutions of the Barkley model (1) are shown for different excitabilities . The spiral wave in Figure 4(a) rotates around a circular core which is small compared to the computational domain. The contour lines of this spiral wave coincide well with a fitted Archimedean spiral (light dashed line; online: green). In comparison, for higher values of , when the core of the spiral wave solution is larger, as depicted in Figure 4(b) (smaller dashed circle; online: red), the contour lines of a spiral wave solution are approximated by an Archimedean spiral (light dashed line; online: green) only sufficiently far away from the core. This effect worsens for larger core radii. The inability of Archimedean spirals to approximate spiral wave solutions near the core can be understood by considering that spiral waves possess locally an approximately constant velocity which is normal to the wave front. However, for uniformly rotating Archimedean spirals the radial velocity is constant. This is particularly problematic close to the origin where the normal direction of the wave is considerably different to the radial direction. Especially for numerical investigations in the large core limit, where numerical domains are not able to include several spiral wave coils, Archimedean spirals will not serve as good approximations within the computational domain.
Note that logarithmic corrections to (30) have been considered (see for example [35]). However, these corrections become negligible in the far field near the boundary, and more importantly for our purpose here as boundary conditions, their derivatives with respect to the will be small at the boundary.

Based on the assumption of a locally normal velocity for spiral waves, it was suggested in [44, 45, 29] to approximate spiral waves by involutes of a circle. In this case contour lines of the spiral wave are given by

(31)

allowing us to simplify . The parameter denotes the radius of the circle that is revolved by the involute. The involute is only defined for . The radius is somewhat arbitrary and we have to choose a proper definition that suits our application. We will later explain how to choose an appropriate . Note that for we have , i.e. the two approximations of Archimedean spiral and involute of a circle coincide in the far field. In Figure 4(b) we show how the approximation of an involute of a circle compares to the one by an Archimedean spiral.

(a)
(b)
Figure 4: Spiral wave solution calculated by solving the stationary frozen system (28) on a circular domain of radius . (a) Spiral solution for with a superimposed Archimedean spiral (dashed line; online: green) with . (b) Spiral solution for with superimposed Archimedean spiral (light dashed line; online: green) with and superimposed involute of a circle (dashed black line) with (dashed circle; online: magenta). The smaller circle (solid line; online: white) with radius is the trace of the tip of the spiral wave defined in (3).

In [33] the two geometric approximations were compared in their ability to approximate experimental data of the Belousov-Zhabotinsky reaction. As the experimentally obtained spiral waves were rotating around small cores, both approximations showed equally good agreement. It is clear that both approximations fail close to the spiral tip [50], however, for our purpose, as illustrated in Figure 4(b), involutes of circles are better suited, especially for the large core limit.

We can use the approximations of an Archimedean spiral or of an involute of a circle, to formulate boundary conditions. Under these approximations we can use contour lines to express spiral wave solutions as

(32)

where is given by equation (30) or (31), respectively. Differentiating (32) leads to

(33)

where is a coefficient depending on which spiral approximation has been chosen. For the Archimedean spiral we find

(34)

Note that is constant for Archimedean spirals. For the involute of a circle we find

(35)

When evaluated at the boundary of the computational domain, we coin this type of boundary condition spiral boundary condition (SBC). SBCs have the advantage that derivatives of a spiral wave solution on the boundary can be expressed entirely by known values of from inside the domain and from the boundary. Note that the geometric approximations can still be used to formulate SBCs, even, if close to the tip, they are not accurate.

The spiral boundary conditions (33) are formulated within the coordinate system which is the coordinate system with origin at the centre of the spiral wave . This does, in general, not coincide with the polar grid of the computational domain. Assume computations are performed on a circular domain of radius in a coordinate system with centre .

To complicate things, the centre of the spiral and the associated coordinate system typically shift during the process of freezing. We therefore need to perform a coordinate transformation relating the two coordinate systems, and , in order to express the spiral boundary conditions (33) in terms of the coordinate system of the computational domain . Using elementary trigonometric relations (see Figure 5) we can write

(36)
(37)

We can now formulate SBCs on the actual computational domain by inserting the transformation into (33), to obtain

(38)

with

(39)

where is given by (34) for Archimedean spirals and by (35) for involutes of circles.


Figure 5: Coordinates of a point in two different polar coordinate systems and .

4.2.1 Determination of the parameters of the spiral boundary condition

The SBC requires knowledge of the location of the centre of the spiral wave and its wavelength in the case of the Archimedean spiral (34), and knowledge of the parameter for the case of the involute of a circle (35). In the following we will explain how these parameters can be determined.

We start with the Archimedean spiral. The freezing method automatically determines the instantaneous centre of rotation of the spiral wave solution via (12). We may therefore determine on the fly during the freezing procedure as

(40)

which becomes exact once the group parameters have become constant. Note that for rigidly rotating spirals unless we are at the bifurcation from rigidly rotating spirals to retracting fingers. Approaching this limiting case, becomes larger, reflecting the increase of the core radius which diverges at criticality. If the tip were fixed at the origin, the core radius could be determined exactly by using only the group parameters. In this case the rotation frequency and the core radius are exactly inversely proportional to each other. In the large core limit, an offset of the spiral tip from the origin of the computational domain, however, becomes negligible as demonstrated later in Figure 14. See [16] where a phase condition is used which pins the tip at the centre of the domain.

The parameter can be estimated on the fly as well during the freezing procedure. For small core spirals we can do so by employing the dispersion relation of travelling wave trains [48]. Before the application of the full D-freezing procedure, the dispersion relation for the velocity of a travelling wave train is determined for each fixed . This is done by applying the freezing method to the Barkley model (1) in a D periodic domain with length . Assuming that the 1D velocity is a good approximation for the normal velocity of the far field spiral wave coils, this velocity should equal the asymptotic far field velocity of a spiral wave , which is consistent with our geometric approximation of Archimedean spirals. The rotation frequency is determined during the freezing procedure and given simply by . The wavelength of the spiral wave is then obtained as the solution of . For sufficiently large core radii, when the spiral wave coils do not interact with each other, we may replace by the 1D velocity of an isolated pulse . The wavelength can then be estimated by . The velocity can also be determined using a D freezing procedure, or via direct simulations of the D version of the Barkley model (1), where the box length is chosen large compared to the decay length of the inhibitor.

For large core spirals, a simpler method can be used to determine . In this limit we can employ the approximation which is independent of or the dispersion relation , which we would need to determine in advance. In Section 5.1 we will present numerical results corroborating these approximations.

The SBC involving the approximation of involutes of a circle requires that we specify the radius of the circle around which the tip revolves. We found that the involute of a circle with radius is a better approximation of contour lines of the actual spiral wave solution than the involute of a circle with the somewhat arbitrary radius , in particular in weakly excitable media with large core radii.

4.2.2 Applicability of the spiral boundary conditions

In the following we investigate the applicability of the spiral boundary conditions (38) and (39). There are two possible problems that can occur: Firstly, the case when there are grid points on the boundary at which but , violating (38). Secondly, the expression for in (39) can possess singularities. We will derive conditions on the location of the spiral centre to ensure that the spiral boundary condition is applicable and both possible problems are avoided. For simplicity, we will use the example of the involute of a circle.

At first we investigate, when there are points on the boundary such that but . This is equivalent to asking whether there exist points on the boundary at which the involute is tangent to the circular boundary of the computational domain. The involute of a circle with radius can be described by

A parametric equation for the circular boundary reads as

Without loss of generality we may set . A necessary and sufficient condition for a tangency at the boundary is given by the solution of the two equations and , which is found to be . Hence, SBCs are violated whenever .

In a next step we investigate the conditions for singularities of . Zeros of the denominator can be found for . Using the trigonometric identity and the coordinate transformation (36), (37) we find , implying as before . We note here without derivation, that for Archimedean spirals one finds [24].

In the small core limit, when the core radius is small compared to the computational domain, the centre of the core can be placed close to the origin of the computational domain, assuring . If the core becomes larger upon increasing , the centre of the spiral wave core will move outside the computational domain, if the spiral wave tip remains resolved within the domain. In this case the conditions can be satisfied, by placing the tip of the spiral wave solution close to the centre of the computational domain, which implies . Recalling our definition , we have

(41)

where we define to be the velocity of the spiral tip tangential to the circle with radius . Because of the positive curvature of the wavefront close to the core, this velocity is smaller than the velocity of the spiral coils in the far field with vanishing curvature. Hence, if the tip is placed close to the centre of the computational domain, is satisfied, assuring that the SBC is neither singular nor violated.

In the large core limit, however, several problems arise, which cannot be resolved by an appropriate positioning of the spiral wave solutions within the computational domain. For the approximation with involutes of circles different problems may arise which need to be addressed. First, for large core spirals the involute circle (or the spiral core) intersects the boundary of the computational domain, and therefore there are parts of the boundary for which , and the SBC (38) with (39) and (35) is not defined anymore. However, in this case we can formulate mixed boundary conditions. We keep SBCs on the part of the boundary where , and we impose the NBC for the remainder of the boundary where , i.e.

(42)

where is given by (39) using (34) for the Archimedean spiral or (35) for the involute of a circle. This is justified since, per definition, the core region of a rigidly rotating spiral is the part of the domain which is never excited by the spiral wave, with and . For the part of the domain which does not lie within the core region, but within the circle of radius , we may also set since in the large core limit, the fields decay fast enough.

Approaching the critical point with , we inevitably reach the point when is larger than the domain size . At this point none of the geometric approximations we discussed is valid anymore. In this case, NBCs, formulated in Cartesian coordinates, prove to be a good approximation, as the curvature of the spiral wave solution becomes negligible and the spiral wave appears to have the shape of a travelling finger.

Analogously to (38) we can formulate spiral boundary conditions for Cartesian grids . We perform a transformation between -coordinates of the spiral system and -coordinates of the computational grid

(43)
(44)

We find

(45)

for boundaries parallel to the -axis, and

(46)

for boundaries parallel to the -axis. Here is given by

(47)

where is given again by (34) for the Archimedean spiral and by (35) for the involute of a circle. The variables and are expressed as functions of using the coordinate transformation (43).

Spiral boundary conditions in Cartesian coordinates can only be applied for large core spirals when the core radius is considerably larger than the box length of the rectangular computational domain. This is due to the fact that on rectangularly shaped domains covering at least one revolution of the spiral arm, there are always points on the boundary at which the spiral arm is tangent to the boundary, violating the SBCs (45) or (46). Therefore we have to restrict the method to large core spirals where only a finger-like part of the spiral is resolved within the computational domain. The orientation of the domain can be chosen such that only one boundary intersects with the spiral wave arm. In this case we can set-up boundary conditions, by imposing an SBC on that boundary and Neumann boundary conditions on the remaining boundaries where activator and inhibitor are assumed to be approximately zero. Without loss of generality we choose to be the boundary which intersects with the spiral. The mixed SBC-NBC boundaries are then written as

(48)

In Figure 6 a contour plot of the activator in Cartesian coordinates is shown, where we used boundary conditions (48). Figure 6 shows an overlay of this solution with a similar solution where computations were performed with NBCs for all boundaries. We show the contour lines of and respectively, and an inset zooming into the area close to the boundary. For the solution obtained by using SBCs, the boundary appears transparent, whereas the solution obtained by using NBCs for all boundaries exhibits a kink near the boundary. The difference between the two solutions is confined near the boundary, and is absent in the interior.

Figure 6: Spiral wave solution in the large core limit obtained by solving the stationary frozen system (27) in Cartesian coordinates. The excitability parameter is , and the spatial discretization is . (a) Contour plot of activator . The inset is added to show the correct aspect ratio. (b) Contour lines (online: black, green) and (online: blue, red) for SBC (involute of a circle) and NBC respectively (inset: zoom into dashed box close to the boundary at ).

4.3 Oscillation-free implementation of boundary conditions

In the middle row of Figure 2 we show results of the stationary problem (28) in polar coordinates using the spiral wave boundary condition (38) for the involute of a circle with (39) and (35). Comparing Figure 2 and Figure 2, we see that the spiral boundary condition respects the shape of the solution at the boundary and does not include spurious kinks. However, as can be seen in Figure 2, the spiral boundary conditions as described in this section, are not able to control the oscillations near the boundary.

In this section we present an explanation for these oscillations, and suggest a simple way to eliminate them. To understand what causes spurious oscillations for NBC and SBC, we investigate their numerical implementation in more detail. Without loss of generality we restrict the discussion here to polar coordinates.

For both boundary conditions, NBCs and SBCs, the equation for the inhibitor involves only first derivatives. This implies, that at the radial boundary , the inhibitor is computed only from values of and on the boundary from the previous time step. Whereas the discrete Laplacian, present in the equation for the activator , couples values at the boundary to those in the interior, i.e. , the boundary values of the inhibitor are decoupled from the interior, and do not receive the outward flowing information from the interior. It is this decoupling of the boundary from the interior for the non-diffusive inhibitor which causes the spatial oscillations near the boundary. (We recall that there are no oscillations for the diffusive activator , and also no oscillations when diffusion is added to the equation for the inhibitor.)

For both cases, NBCs and SBCs, we propose a simple method to overcome this decoupling. For the activator we use NBCs or SBCs as discussed. However, instead of invoking this boundary condition for the inhibitor as well, we evaluate the derivative of the inhibitor at the boundary by a one-sided second-order discretization according to

(49)

For Cartesian coordinates we use equivalent expressions.
Considering the original unfrozen Barkley system (1), there are two reasons why this simple trick works. Firstly, the boundary condition is consistent with the fact that information flows outwards when studying spiral waves, and therefore outer grid points are influenced by inner grid points only. Secondly, the inhibitor is “slaved” to the activator in the sense that if the activator is known, the linear equation for the inhibitor can be integrated explicitly. The boundary conditions of the activator are therefore (approximately) inherited by the inhibitor.

In Figure 2 (bottom row) we show how the implementation of the spiral wave boundary condition for the activator and the one-sided derivative (49) for the inhibitor in the stationary problem (28) in polar coordinates suppresses the spatial oscillations and avoids changes in shape and amplitude (see also Figure 3). In Figure 7 we show how the implementation of Neumann boundary conditions for the activator and the one-sided derivative (49) for the inhibitor controls the spatial oscillations and the spiral wave shape in the stationary problem (27) in Cartesian coordinates in the weakly excitable regime.

During extensive numerical tests [24] we found that in highly excitable media only spiral boundary conditions (38) in conjunction with the one-sided boundary condition (49) give accurate results for the stationary and the time-dependent freezing method without spurious oscillations and without unnatural deformations of the spiral shape. For the large core limit, in weakly excitable media, Neumann boundary conditions in Cartesian coordinates are the optimal choice to avoid numerical oscillations and to properly resolve the shape of the spiral wave solution. Note that in the large core limit NBCs resemble SBCs. This is due to the fact that the overall curvature of the spiral wave, captured in the computational domain, is close to zero. In a way NBCs can be viewed as a simpler and more convenient formulation of SBCs in the large core limit, which do not involve the core radius and the wavelength of the spiral wave.
In the following we will be using SBCs in polar coordinates for highly excitable media, and NBCs in Cartesian coordinates in the large core limit. We will also employ the one-sided boundary condition (49) in all numerical simulations to avoid oscillations.

Figure 7: Inhibitor in the weakly excitable case with , calculated from the stationary frozen system (27) in Cartesian coordinates with NBCs for activator and inhibitor (top panels) and NBCs for activator only and one-sided derivative (49) for the inhibitor (bottom panels). (b) and (d) are close-ups of the respective solutions depicted in (a) and (c), zooming in on the boundary at .

5 Numerical investigation of the large core limit of spiral waves

In this Section we provide numerical results on the wavelength , the core radius and the rotation frequency of spiral waves in the large core limit. The rotation frequency and the core radius are computed from the group parameters . We recall that and can be calculated using (13) with the centre of the spiral given by (12) and the location of the tip defined as the intersection of the contour lines (3). We approach numerically criticality where the core radius develops a singularity and the frequency approaches zero. Before we present these results, we collect and discuss some practical aspects of the different methods we introduced to study spiral waves in the large core limit.

We solve the stationary freezing problem (27), increasing . As initial guess for the Newton-Raphson method we use the frozen solutions of the previous value for . We use a discretization of for Cartesian simulations, unless stated otherwise. We have checked the quadratic convergence of the error in the calculation of the core radius and the rotation frequency with respect to the spatial discretization, and found that at this resolution the values have sufficiently converged.

For polar coordinates we use a discretization and . For a computational domain with this corresponds roughly to an equivalent discretization of and near the boundary. This immediately alludes to practical limitations of polar coordinates for large computational domains; wave profiles away from the centre of the computational domain are not properly resolved, and has to be sufficiently small, unless one employs a computationally expensive fine angular discretization. At this point it is important to repeat, that in highly excitable media with small core radii at sufficiently small values of , polar coordinates work very well. This is due to two effects: First, transverse wave profiles are wider for smaller values of , and second, the smaller wavelength in this case causes a cross section of the spiral wave at the boundary which is wider than its transverse profile.

In Figure 8 we show results for solving the stationary frozen system (28) in polar coordinates. We used two types of boundary conditions, i.e. NBCs (red ) and SBCs based on the approximation by involutes (blue ). We further included results from direct simulations of the full Barkley model (1) (black ). For a computational domain of size direct simulations are limited to excitabilities of with corresponding radii . One can, in principle, determine the core radius for parameter values larger than , however, only with the additional computational cost of increasing the computational domain. There will be a value of where this is not possible anymore with given computational power. We stress that this value of is way below what we call the large core limit. This can be seen in Figure 9 (left panel) where we show contour plots of the activator at different values of corresponding to data points in Figure 8 for SBCs. The white circular lines represent the trace of the tip as defined in (3) and are calculated by solving equations (17) for the group variables. Figure 9(c) indicates that it becomes inherently difficult to calculate spirals and their characteristic parameters and for large cores by direct simulations of the full Barkley system (1). It is this observation which makes the freezing method so attractive for large core spirals.
The insets in Figure 8 clearly show that SBCs outperform NBCs in the highly excitable regime, and allow for the determination of frozen solutions and the values of the core radius and the rotation frequency for a greater range of parameter values. The breakdown of NBCs is linked to the finite size of the computational domain as illustrated in Figure 10. We can see that the finger is not properly oriented in the finite domain and that therefore NBCs are not a suitable choice, bending the solution into an unnatural direction. This does not happen for the shape-preserving SBCs. The negative effect of the (inaccurate) NBCs to deform contour lines of the solution is proportionally larger for moderate computational domains than for larger ones (which makes polar coordinates more sensitive to these effects then Cartesian coordinates). In principle, one may extend the range of validity in -parameter space for NBCs by considering larger and larger computational domains; however, this would quickly become computationally unfeasible.

Figure 8: Core radius and rotation frequency of spiral waves as functions of computed in polar coordinates. Here ‘’ denotes values obtained by a direct simulation of the Barkley model (1), and ‘’ and ‘’ represent results from solving the stationary frozen system (28) with NBCs and SBCs respectively. A computational domain with was used.
(a) (b) (c) (d) (e) (f)
Figure 9: Contour plots of the activator solution of the stationary frozen system in polar coordinates with SBCs based on approximations by involutes (left) and in Cartesian coordinates with NBCs (right) for increasing values of excitability . The white circular lines indicate the trace of the tip.
(a) SBC
(b) NBC
Figure 10: Contour plot of the activator found as the solution of the stationary frozen system (28) at with with SBCs (left) and NBCs (right), respectively.

To approach the large core limit, we employ a Cartesian coordinate system for values of . In Figure 11 we show the rotation frequency and the core radius of simulations where we approach the critical point at . Here we use and . In Figure 9 (right panel) we show contour plots of the activator solution of the frozen stationary system (27) at different values of corresponding to data points in Figure 8 for NBCs. The exact value of depends on the discretization and also on the size of the computational domain. See Section 5.2 for a discussion on this issue. Compare the range of attainable in polar coordinates and Cartesian coordinates (cf. Figure 8 and inset of Figure 11). Contrary to the results with polar coordinates, in Cartesian coordinates NBCs are applicable for a larger range of parameter values than SBCs. As discussed in Section 4.2.2, SBCs are not applicable in the large core limit, when the spiral wave appears as a finger in computational domains of finite size. The part of the spiral wave solution, starting at the tip, which cannot be approximated by an Archimedean spiral or an involute of a circle, increases the closer one is to criticality. As a proxy for the extent of this region we depict in Figure 12, how grows with . For values , corresponding to core radii , we have , and the freezing method using SBCs is not applicable anymore. NBCs, on the other hand, become a better approximation the closer we are to the critical point, where the spiral wave solution approaches zero curvature – provided that the finger is appropriately oriented within the computational domain to assure a perpendicular intersection with the boundary.

To study the behaviour of the wavelength, the rotation frequency and the core radius of spiral waves in the large core limit, we therefore use from now on Cartesian coordinates and Neumann boundary conditions. In Figure 13 we show contour plots of the activator and the inhibitor close to criticality, illustrating the appropriateness of Neumann boundary conditions in Cartesian coordinates in this limit.

Figure 11: Core radius and rotation frequency of spiral waves as functions of with computed in Cartesian coordinates. Here ‘’ denotes values obtained by direct simulation of the Barkley model (1), and ‘’ and ‘’ represent results from solving the stationary frozen system (27) with NBCs and SBCs respectively. A computational domain with and was used.
Figure 12: Core radius and radius of the circle of the involute as a function of . The critical value of the excitability parameter is .
(a) Activator
(b) Inhibitor
Figure 13: Contour plot of a spiral wave solution of system (27) at , close to the critical excitability . The (a) activator and (b) inhibitor are shown with the inset depicting them with the correct aspect ratio. A computational domain with and was used.

5.1 Wavelength

The application of the spiral boundary conditions using Archimedean spirals requires the knowledge of the wavelength . In the small core limit the wavelength can be determined as the solution of the implicit equation

(50)

where is the velocity of a 1D wave train with wavelength (see Section 4.2.1). In the large core limit, where the inhibitor decays sufficiently quickly and spiral wave coils do not interact, we may further simplify to , where is the velocity of an isolated 1D pulse. In Figure 14 we show a comparison of these expressions with numerical results from a direct simulation of the full Barkley model (1). The velocities and are determined by freezing pulses in the corresponding 1D-model with box length using the same discretization as in two dimensions. We see that the wavelength determined by (50) is a reasonably good approximation of the true wavelength even for small core radii. For larger radii the two methods to determine converge.
In the large core limit we can deduce a simpler approximation for the wavelength which does not require the independent determination of the 1D velocities. One can define two approximate temporal periods for rigidly rotating spiral waves with wavelength and curvature . First, the temporal period of a spiral wave with velocity given to first approximation by , and second, the time which measures the time of one revolution of a spiral wave tip around a circle with radius chosen such that the normal velocity of the spiral tip is tangential to that circle. Equating these two temporal periods leads to the kinematic relation [26]

(51)

In the large core limit, we expect , and