Unstable Spiral Waves and Local Euclidean Symmetry in a Model of Cardiac Tissue
Abstract
This paper investigates the properties of unstable single-spiral wave solutions arising in the Karma model of two-dimensional cardiac tissue. In particular, we discuss how such solutions can be computed numerically on domains of arbitrary shape and study how their stability, rotational frequency, and spatial drift depend on the size of the domain as well as the position of the spiral core with respect to the boundaries. We also discuss how the breaking of local Euclidean symmetry due to finite size effects as well as the spatial discretization of the model is reflected in the structure and dynamics of spiral waves. This analysis allows identification of a self-sustaining process responsible for maintaining the state of spiral chaos featuring multiple interacting spirals.
pacs:
02.20.-a, 05.45.-a, 05.45.Jn, 47.27.ed, 47.52.+j, 83.60.WcUnstable solutions of reaction-diffusion models often inherit the Euclidean symmetries of the underlying evolution equations despite the presence of boundary conditions, or spatial discretization, that break those symmetries. A good example is provided by spiral traveling waves, which arise in models of excitable systems possessing both translational and rotational invariance. In particular, unstable spiral waves are believed to play a major role in initiating and sustaining cardiac arrhythmias such as atrial or ventricular fibrillation. While stable spiral solutions can be readily computed in any geometry, computing unstable spirals can be rather non-trivial, especially when the cores of the spirals drift. We show how unstable spiral waves can be computed on domains of arbitrary shape without making any simplifying assumptions and investigate how the structure and dynamics of these waves depends on the domain size, the position of the spiral core, as well as the relevant parameters of cellular kinetics.
I Introduction
Spatially coherent contractions that occur in the heart are essential for pumping blood through the circulatory system. The atria, followed by the ventricles, relax, filling with blood, and then contract in a spatially coherent way, pushing the blood out. This coordinated relaxation and contraction, referred to as normal rhythm, can be disrupted by various mechanisms, leading to a range of arrhythmic behaviors, in which the efficiency of the heart as a pump is reduced, sometimes drastically. Fibrillation is the most complex regime characterized by “turbulent” dynamics featuring multiple interacting unstable spirals where both spatial and temporal coherence of the contractile dynamics is destroyed. Ventricular fibrillation Cherry and Fenton (2008); Gray, Pertsov, and Jalife (1998) is particularly dangerous and, if not treated within minutes, causes death.
Numerous models of cardiac tissue of varying complexity exist. The majority of so-called monodomain models fall in the general class of reaction-diffusion systems Clayton et al. (2011); Courtemanche (1996); Barkley (1991); Karma (1994); Bueno-Orovio, Cherry, and Fenton (2008); Simitev and Biktashev (2006); Noble (1962); Cherry et al. (2007). Besides cardiac tissue, reaction-diffusion systems are also used to model diverse phenomena such as chemical reactions Muller, Plesser, and Hess (1985); Keener and Tyson (1986), bacterial chemotaxis Siegert and Weijer (1991); Vasiev, Hogeweg, and Panfilov (1994), and disease propagation Murray, Stanley, and Brown (1986). All of these systems display turbulent solutions dominated by multiple interacting spiral waves, though the terminology varies: such solutions are referred to as defect-mediated turbulence Ouyang and Swinney (1991), spiral chaos Hagberg and Meron (1994), spiral breakup Panfilov and Hogeweg (1993), or spiral defect chaos Morris et al. (1993). This lack of consistency reflects the mostly empirical approach to the study of multi-spiral patterns and our lack of fundamental understanding of their dynamics.
Previous attempts to build a dynamical description of fibrillation (and more generally, spiral chaos) were based mostly on the intuition gleaned from the studies of stable solitary spirals. There is a vast literature on the subject, so we will only mention studies most relevant to the present investigation. In particular, Barkley et al. Barkley, Kness, and Tuckerman (1990) showed that even very simple reaction-diffusion models can produce qualitatively different types of spirals. The spiral tip can move either in a circular trajectory or in a periodic or quasi-periodic fashion. By using a reference frame rotating with the spiral, the tip dynamics can be simplified: in the first case the tip becomes stationary while in the second it executes a periodic motion. Respective spirals are referred to as pinned or meandering. Our use of the term “pinned” to describe non-drifting spiral waves deviates from the common usage – in which a spiral wave is attached to macroscopic heterogeneities of the tissue or medium – only in scale. Since this distinction is superficial, we will use this term to describe spirals pinned to microscopic features as well. Barkley Barkley (1992) subsequently showed that Newton-Krylov method can be used to efficiently compute pinned spiral waves in a rotating reference frame in which they become stationary (i.e., are described by relative equilibria). He also computed their leading eigenvalues using Arnoldi method Goldhirsch, Orszag, and Maulik (1987) and verified that transition to meandering spirals (relative periodic orbits) is described by a Hopf bifurcation. The same approach has later been applied for computing unstable spiral waves in a model of cardiac tissue by Allexandre and Otani Allexandre and Otani (2004).
Barkley Barkley (1994) and Barkley and Kevrekidis Barkley and Kevrekidis (1994) showed that the tip dynamics for both pinned and meandering spirals can be reproduced using a system of five weakly nonlinear ordinary differential equations (ODE) derived by assuming the dynamics are equivariant under a Euclidean symmetry group. Fiedler et al. Fiedler et al. (1996); Fiedler and Turaev (1998) showed that, more generally, equivariance with respect to noncompact, finite-dimensional Lie groups (such as the Euclidean group ) allows description of the dynamics near relative equilibria (such as rigidly rotating spirals) in terms of a skew-product flow, where the motion transverse to the group manifold is decoupled from the motion on the group manifold. Locally, the group manifold represents all symmetry transformations of a particular solution. For instance, in two dimensions, the skew decomposition separates the evolution of the shape of the spiral from the changes in the position or phase of the spiral. Sandstede et al. performed a center manifold reduction of the dynamics near relative equilibria Sandstede, Scheel, and Wulff (1997) and near relative periodic orbits Sandstede, Scheel, and Wulff (1999), formalizing and extending the reduced description of Barkley and Kevrekidis Barkley (1994); Barkley and Kevrekidis (1994) to physical spaces of arbitrary dimensionality. An in-depth discussion of the role of symmetries in dynamics is available in, e.g., Chossat and Lauterbach Chossat and Lauterbach (2000).
Beyn and Thümmler Beyn and Thümmler (2004) developed a numerical method for computing the dynamics near rotating relative equilibria on unbounded domains which used the skew-product representation of the dynamics to eliminate or “freeze” the dynamics along the group manifold. The freezing approach was later used by Beyn and Lorentz Beyn and Lorenz (2008) to numerically compute the entire stability spectra for pinned spiral waves. They also found good agreement between the numerically computed eigenvectors associated with marginal eigenvalues and the Goldstone modes associated with infinitesimal translations and rotations of the spiral wave. The same approach was later used by Hermann and Gottwald Hermann and Gottwald (2010) to investigate the dynamics of spiral waves in the large-core limit and by Foulkes and Biktashev Foulkes and Biktashev (2010) to characterize drift and meandering of spiral waves.
In the case of fibrillation, individual spirals possess rotational frequencies in excess of the normal rhythm pacing Fenton, Cherry, and Glass (2008) and as a result are typically strongly unstable, often encountering refractory regions of tissue and breaking up within a few rotations Bär and Eiswirth (1993); Fenton et al. (2002); Bernus, Verschelde, and Panfilov (2003). However, at present our understanding of the properties and dynamics of unstable spirals, especially in the context of cardiac dynamics, is limited. This study investigates unstable single-spiral wave solutions in a simple model of cardiac tissue. The key objectives are to understand how their properties are affected by the size of the computational domain, the proximity of the boundaries, and by the spatial discreteness, which is an essential feature of cardiac models reflecting the finite size of cardiac cells – cardiomyocytes.
The paper is organized as follows. We begin in Sect. II with an overview of reaction-diffusion systems and, in particular, the Karma model of cardiac tissue used in the remainder of this study. The numerical method used to find exact unstable non-chaotic solutions is described in Sect. III. The results of numerical experiments investigating the properties of unstable single-spiral solutions of the Karma model and the emergence of local Euclidean symmetries are presented in Sect. IV. Sect. V summarizes these results and discusses their implications in the context of fibrillation (or, more generally, spiral chaos).
Ii The model
Although our focus here is primarily on the dynamics of cardiac tissue, much of the subsequent discussion applies in equal measure to a broader class of excitable systems Biktashev (2007) and, even more generally, to reaction-diffusion systems that support wave propagation. Reaction-diffusion equations were originally introduced to describe pattern formation arising in systems involving reacting and diffusing chemicals. If one forms a vector field from the concentrations of interacting chemicals, , their evolution can be described by the system of partial differential equations (PDE)
(1) |
where is a diagonal matrix of molecular diffusion coefficients, and the two terms on the right-hand-side describe, respectively, diffusion and reactions.
Similar equations, however, have also been used to describe physical and biological systems that do not necessarily involve molecular diffusion. In particular, monodomain models of cardiac tissue dynamics have the same form, although the interpretation of both the variables and the terms on the right-hand-side is different. To distinguish the physical nature of different variables, we will introduce a different notation, , where is an (electric) transmembrane potential and is a set of gating variables.
Due to the variation between species and even between the different regions of the same heart (e.g., atria, ventricles, Purkinje fibers), there is not a universal model of cardiac dynamics akin to the Navier-Stokes equations for fluid flow. While the equation describing the transmembrane potential, at least in the monodomain formulation, is universal, the models of cellular kinetics vary widely in complexity and even the number of gating variables Fenton and Cherry (2008). Physically more realistic bidomain models Sepulveda, Roth, and Jr. (1989) of two- and three-dimensional tissue include an additional nonlocal constraint equation (for the extracellular potential), which makes them significantly more expensive computationally. However, their dynamics are qualitatively similar to those of the monodomain models, and when the anisotropy of the intracellular and extracellular media is the same, bidomain models can be reduced to the monodomain formulation Roth (2008). Hence, our focus here will be entirely on the monodomain model of cardiac tissue dynamics.
One of the simplest models, which reproduces dynamics qualitatively similar to fibrillation, is due to Karma Karma (1994). The Karma model was formulated as a minimal model of cardiac dynamics exhibiting alternans – the instability that is believed to be responsible for initiating and sustaining fibrillation in the heart Karma (1993, 1994). We use the following nondimensionalized form of the “reaction” terms:
(2) |
where there is only one gating variable, so that .
Several modifications were made to the original model in the interest of creating a well-conditioned dynamical system. All instances of the Heaviside step function were replaced, following Allexandre and Otani Allexandre and Otani (2004), with a smooth function . This substitution removes the singularities in the partial derivatives while recovering the original step function in the limit . An additional term was introduced to smooth out an unphysical singularity in the variable which forms at the spiral core. In the absence of this term, numerical solutions of the original model become extremely sensitive to the spatial discretization. Finally, the original model assumes no diffusion for the gating variable. In the modified system, we set with , to further regularize the dynamics in the core region. Although diffusion is usually ignored in the equations describing the dynamics of gating variables, all relevant ions and even most secondary messengers such as IP, cAMP and cGMP can pass through the gap junctions between cells Bevans et al. (1998); Garcia-Dorado, Rodriguez-Sinovas, and Ruiz-Meana (2004), so a small nonzero value of is justified from the physiological perspective. These modifications have a very subtle effect on the cellular dynamics. In particular, the -nullcline (Fig. 1a) is smoothed into a sharply varying, but -continuous function (Fig. 1b). The position and the stability of the equilibria remain essentially unchanged for sufficiently large, however, significant deviation in the position of the equilibria is found as is decreased to .
We kept the parameter values used in the original work Karma (1994), i.e., , , and . All times and lengths were non-dimensionalized using, respectively, the characteristic time scale ms and the characteristic length scale cm (the size of a typical cardiomyocyte). The original diffusion constant cm/s was non-dimensionalized accordingly, yielding . The ratio of the two diffusion constants was set to . Finally, the width of the smoothed step function was varied in the range .
The dynamics of the modified model with are essentially indistinguishable from those of the original model, so we can rely on the original analysis in choosing the value of the restitution parameter . Karma showed Karma (1994) that the pinned spiral solution is stable at intermediate values of . As is decreased below about 0.8 (0.6 according to our more precise calculations), it becomes unstable and gives rise to a stable meandering spiral via a Hopf bifurcation. If is increased above about 1.0, the pinned spiral also undergoes a Hopf instability, which however does not lead to meander. Instead, one finds temporal and spatial modulation of the spiral wavelength (or action potential duration), which is referred to as alternans in the cardiac literature Guevara et al. (1984); Nolasco and Dahlen (1968). As is increased further, conduction block occurs away from the spiral core, leading to a breakup of the spiral and eventual transition to the state of spiral chaos featuring multiple interacting spirals. In the remainder of this paper we set (which corresponds to ). For this value of spiral chaos persists indefinitely for sufficiently large domains Karma (1994).
We used the physiologically and dynamically relevant no-flux boundary conditions, . These boundary conditions properly describe the vanishing of electric current at the boundaries of cardiac tissue (or the boundaries of its excitable regions). In addition, our investigation Byrne, Marcotte, and Grigoriev (2015) of multi-spiral states characterizing spiral chaos showed that at the edges of domains (or tiles) supporting individual spirals.
ii.1 Spiral Waves and Symmetry
Reaction-diffusion systems (1) with arbitrary local kinetics are equivariant with respect to Euclidean transformations , where for an -dimensional physical space, i.e., . The subgroup reflects the invariance with respect to time translation, but not time-reversal. In an unbounded two-dimensional domain, the corresponding symmetry group includes continuous translations along the two coordinate directions , continuous rotations in the plane, a discrete inversion symmetry in space, and continuous forward-time translations. As a result, reaction-diffusion systems, and the Karma model in particular, possess solutions that respect some or all of these symmetries on spatially unbounded domains.
On unbounded domains and in the absence of discretization these symmetries are exact. As a result there are solutions, known as relative equilibria, for which time-evolution is equivalent to a symmetry group transformation. Due to this special property, relative equilibria can be reduced to an equilibrium solution in a moving (e.g., translating or rotating) reference frame. One such example is a rotating, or spiral, wave which satisfies , where is the angular speed and is the angular derivative about an origin . In this particular case, time evolution is equivalent to a rotation about that origin with angular speed .
Another, more complex, class of spiral wave solutions is described by relative periodic orbits, i.e., solutions which repeat exactly after some time , , up to a symmetry transformation . This transformation can include rotations (for meandering spiral waves) or translations (for drifting spiral waves). Similar to relative equilibria, in an appropriate (rotating or translating) reference frame a relative periodic orbit can be reduced to a time-periodic solution.
Spiral wave solutions described by relative equilibria and relative periodic orbits can only be found on unbounded domains, or domains with periodic boundary conditions, since generic boundary conditions break both translational and rotational symmetries (one exception is solutions which respect rotational symmetry relative to the center of a circular domain). Since the motivation for this study is to better understand multi-spiral states, the ability to compute single-spiral solutions on domains of arbitrary shape and size is essential. Furthermore, cardiac tissue is heterogeneous, with the cellular-level models possessing an intrinsic length scale defined by the size of cardiac cells. This microscopic heterogeneity also breaks Euclidean symmetry of the PDE (1).
It is fairly well understood how stable spiral wave solutions change once boundaries are imposed. Spiral wave interaction with boundaries, heterogeneities, etc. is controlled by the size of its core, which is determined by the spatial extent of the response functions (adjoint eigenfunctions with unit Floquet multipliers) Biktashev, Holden, and Nikolaev (1996); Biktasheva et al. (2009); Biktasheva, Biktashev, and Foulkes (2010). In particular, if the distance between the tip of the spiral wave and the boundary is larger than , the structure and dynamics of the wave remain essentially unchanged compared with the corresponding solution on an unbounded domain. However, because of symmetry breaking, a pinned spiral wave would be formally described by a periodic solution, rather than a relative equilibrium, in the presence of boundaries. Similarly, a drifting or meandering wave generally would not be described by a relative periodic orbit in the presence of boundaries. We can expect the same conclusions to also apply to unstable spiral wave solutions.
These complications limit the usefulness of global symmetry reduction as a method for computing spiral wave solutions of cardiac tissue models or characterizing their properties, especially in the unstable regime. Therefore, we developed a different, more general, procedure for computing unstable spiral wave solutions on bounded domains of arbitrary shape by reducing the symmetries locally. It is described in the next section using square domains as a fairly representative example.
Iii Computation of unstable spiral waves
Cardiac tissue is made up of muscle fibers which are, in turn, made up of cardiomyocytes. This structure breaks, perhaps weakly, the rotational symmetry of the system. Hence, even though the choice of coordinate directions in the original PDE (1) is arbitrary, the cellular structure imposes a natural choice on the coordinate directions. To be specific, we will assume that the () axis is oriented along (transversely to) the fibers.
The modified Karma model (1)-(2) is solved numerically on square domains of side length by using a finite-difference discretization of the Laplacian operator on a two-dimensional grid of size . Since we chose the size of cardiac cells as the characteristic length scale, , by construction the grid spacing is and in non-dimensional units. (In reality the cells are not square, but we will ignore this complication, along with the differences in the diffusion constant along and transversely to the fibers, since these can be accounted for by simply choosing different length scales in the two coordinate directions.) The no-flux boundary conditions were enforced using the well-known “ghost point” method.
To minimize the impact of discretization on the structure and stability of solutions, the Laplacian operator acting on the state at position is approximated using a second-order finite-difference stencil
(3) |
with coefficients , , and . For , this expression is the most isotropic formulation of the nine-point finite-difference stencil on a two-dimensional uniform grid Lindeberg (1990), and is thus expected to most faithfully recover the rotational and translational symmetry of the original PDE.
The resulting spatially discretized equations are solved numerically using the classical fourth-order Runge-Kutta method, with nondimensional time step . This time step is sufficiently small to accurately resolve the fast time scales in the evolution of the field, and found to produce solutions accurate to one part in in the -norm for integration times of order (or 100 ms, in dimensional units) with the expected convergence, which is sufficient for computing unstable solutions on time scales relevant for the dynamics of this system.
Direct numerical simulation (DNS) can only be used to compute stable (or marginally stable) solutions. Unstable spiral waves described by relative equilibria have been computed with the help of Newton-Krylov method by using rotating reference frames to convert them to equilibria Barkley (1992); Allexandre and Otani (2004). However, this approach cannot be extended to finite noncircular domains on which spiral waves are not described by relative equilibria. Recently a new generation of Newton-Krylov methods was developed Viswanath (2007) that enables efficient computation of many common types of unstable solutions of spatially extended systems, including equilibria and periodic orbits, as well as their relative counterparts in the presence of global continuous symmetries. Specialized versions of these methods applicable to systems with local continuous symmetries are discussed next.
iii.1 Newton-Krylov Solver
We will illustrate the Newton-Krylov method using relative periodic orbits on an unbounded domain as an example. Consider a drifting spiral wave solution which exactly recurs after a period , but shifted by relative to the initial condition (generalization to the case of meandering spiral waves involving rotation is straightforward). Using the evolution operator and the operator of spatial translations , this can be written as . This solution can be specified completely by a real-valued vector , where the state is composed of the values of the fields and on the computational grid. More generally, a point that lies on a relative periodic orbit corresponds to a root of the vector function ,
(4) |
which measures the difference between the “initial” state and the transformed “final” state (here, the more general group transformation takes the place of the translation operator ). Absolute (non-relative) periodic solutions are merely a special case: they satisfy (4) with , (i.e., ). Similarly, relative equilibria may be treated as special cases of periodic orbits.
Taylor expansion of (4) shows that the leading-order approximation for the correction to the system state is given by the solution of
(5) |
where . However (5) does not determine the corrections specifying the group transformation and the periodicity, (e.g., the and of our drifting spiral). This is a generic feature of systems with continuous symmetries, since there is a continuum of solutions satisfying (4). Additional constraints are needed to uniquely specify the solution. Specifically, we will require that the correction to the state be transverse to the group manifold. The constraints can be written using the generators of translation along , , and , which are given by the respective partial derivatives, yielding the following well-posed linear problem:
(6) |
where , and
(7) |
where denotes the final state, and is the row vector of partial derivatives with respect to the independent coordinates which define the local tangents of the group manifold, e.g., for the drifting spiral.
The linear system (6) effectively implements the skew-product representation of the dynamics, where the correction is split into components transverse to, and along, the group manifold. Many authors Knoll and Keyes (2004) drop the dependence on the final state altogether, instead substituting . We are unaware of any benefit in doing so and have opted not to, as is explicitly available: . It should be noted that the transversality condition with respect to rotations is automatically satisfied, since for relative periodic orbits rotations are equivalent to a combination of spatial and temporal translations, obviating the need for an explicit constraint on the “phase” of the solution Barkley (1992). For example, in the special case of relative equilibria we have .
Typical discretizations of (2) involve variables, making it very expensive to even compute the elements of the matrix , let alone solve the system (6) directly. Instead, the solution can be found approximately using a truncated spectral representation of the matrix, which is constructed using an orthogonal basis produced by Arnoldi iteration Goldhirsch, Orszag, and Maulik (1987),
(8) |
where the action of is approximated using finite-differencing with some and denotes the normalized projection of onto the orthogonal complement of the linear space spanned by .
Projecting (6) into the Krylov subspace spanned by we obtain a -dimensional linear system
(9) |
where is a real-valued upper-Hessenberg matrix generated explicitly during the Arnoldi iteration, , , , , and is the projection operator. The vector defines the shift and period corrections in our drifting spiral example. Due to strong contraction associated with the Laplacian term in (1), accurate approximation of the original system (6) can usually be obtained with of order a few tens, very small compared with the full dimensionality of the discretized system, yielding the solution .
Newton-Krylov iterations converge provided a good initial condition is available. In order to improve robustness of the procedure for initial conditions that are further away from the solution, a line search was implemented for the cases when the reduction in residual norm was smaller than the expected improvement from linearization, . A scaling factor for the Newton step magnitude was computed by minimizing the 2-norm of the residual , yielding a more robust sequence of iterations . The iteration is considered converged to the solution when the 2-norm of the residual function is minimized below a predetermined value, (in this study we used ).
For pinned single-spiral solutions on arbitrary bounded domains described by absolute periodic orbits the implementation of the Newton-Krylov method described above is sufficient. In fact, it can be simplified by setting and dropping the condition of transversality with respect to translations, so that . However, Newton-Krylov iterations predictably fail to converge for drifting spirals for which . Since boundary conditions break translational symmetry, the residual (4) cannot be reduced below for some positive constant .
iii.2 Weighted Newton-Krylov Solver
Boundary conditions do not merely break the translational and rotational symmetries. Finite rotations on a bounded domain are not injective: some points are mapped out of the domain and others into (Fig. 2a). A similar situation is encountered for finite translations on a bounded domain (Fig. 2b). Consider, for example, a meandering spiral wave for which on an unbounded domain. On a bounded domain the residual will not vanish (we set outside of to make the residual well-defined). If one places the origin of rotation near the tip of the spiral wave, the residual will decompose into two easily identifiable contributions. Inside (the octagonal overlap region in Fig. 2a), the residual is small and lies near the group manifold,
(10) |
where describes the magnitudes of rotations or shifts accounting for the arbitrary choice of the origin and the frequency dependence on the domain size. Outside the overlap region (the eight triangular regions in Fig. 2a), the residual is large, .
Figure 2b describes the effect of boundaries on drifting spiral waves described by relative periodic orbits for which on an unbounded domain. On a bounded domain we find a decomposition of the residual analogous to the case of meandering waves. In the overlap region (the rectangular region in Fig. 2b), the residual again can be represented in the form (10), where the small shifts and rotations account for the dependence of the drift velocity and rotation frequency on the domain size. Outside of the overlap region, again the residual is large, . Therefore, our objective is to minimize the residual inside the overlap region and suppress it everywhere else.
Generalized relative periodic orbits on bounded domains can be found using a modification of the traditional Newton-Krylov method which introduces an auxiliary weighting (or windowing) of the residual,
(11) |
where is a diagonal matrix with elements . The first diagonal elements correspond to the weights associated with the dynamical field variables and ; they correspond to the values of the windowing function
(12) |
where denotes the center of the domain and determines the diameter of the (circular) “window” in units of (we set and , unless specified otherwise). Specifically, . The remaining diagonal elements correspond to the spatial displacements (if applicable) and the period of the spiral wave and are all set to unity. The effect of windowing on the residual is illustrated for the cases of a meandering and a drifting spiral in Fig. 2c and Fig. 2d, respectively. We should point out that the weighting approach is not limited to square domains and rectangular grids and can be easily applied to domains of any shape with any grids, including unstructured ones.
We should also point out the closely related applications of weighting functions in numerical methods for PDEs such as the phase-field boundary method Bueno-Orovio, Perez-Garcia, and Fenton (2006), domain decomposition Bjorstad and Gropp (2004); Xu (1992), and the damping filter method Teramura and Toh (2014). Alternatively, the weighting may by interpreted as an ad-hoc Jacobi-type preconditioning method for solving linear problems Trefethen and Bau (1997). Even though the use of weighting is not a novel numerical approach per se, to the best of our knowledge, it has never been used either for restoring broken symmetries or for computing unstable traveling wave solutions.
The linear system (11) can be solved in the same way as (6). We found that the use of weighting dramatically improves the robustness of the Newton-Krylov method, not only allowing computation of generalized relative periodic orbits, but also substantially improving convergence speed for (absolute) periodic orbits. In practice, the convergence properties of the weighted Newton-Krylov solver were found to be fairly insensitive to the shape of the function , provided that it vanishes near all the boundaries. In particular, solutions which satisfy (6) can be computed using (11) combined with the relaxation process in which (or ). This idea is similar to the damping filter method Teramura and Toh (2014).
Whether weighting is used or not, the Newton-Krylov solver generates the Floquet multipliers and Floquet modes of the computed solution, i.e, eigenvalues and eigenmodes of the full-space Jacobian ,
(13) |
essentially “for free.” Indeed, the spectrum of the Krylov-subspace Jacobian yields a good approximation to the leading eigenvalues and eigenmodes of , while . Of particular importance are the unstable modes (which correspond to ) and the Goldstone modes (which correspond to ) that characterize the symmetries of the system.
Iv Results and Discussion
In this Section we present unstable spiral wave solutions of the modified Karma model computed using the weighted Newton-Krylov method described in Sect. III and investigate their properties. In particular, the spiral wave solution shown in Fig. 3a was found using continuation of a stable single-spiral solution with an approximately centered tip , achieved by decreasing (which corresponds to increasing the restitution parameter ). It has a period and wavelength , or cm in dimensional units. As Fig. 3b shows, the solution has one Goldstone mode and three unstable modes.
The Goldstone mode presented in Fig. 3f corresponds to the temporal derivative of the solution , as expected for a periodic orbit. The real unstable mode (Fig. 3g) corresponds to an almost rigid shift of the spiral wave in the direction and can be identified with a frustrated translational Goldstone mode (its amplitude varies in space). The stable eigenvalue also corresponds to a frustrated Goldstone mode (Fig. 3e), which can be identified as linear combination of translations in the and directions, , where is the direction of translation. The two complex conjugate unstable modes (Figs. 3c and 3d) correspond to the variation in the width of the excitation wave (i.e., alternation of the action potential duration), which is to be expected for . Since the Goldstone modes associated with spatial translations are frustrated for (continuous translational symmetry is broken), no additional constraints beyond orthogonality with respect to are needed. The residual is minimized for , so this spiral wave is pinned and corresponds to an absolute periodic orbit.
Our results should be contrasted with those obtained by Allexandre and Otani Allexandre and Otani (2004) for the modified version of the 3-variable Fenton-Karma (3V-FK) model Fenton and Karma (1998) (they also replaced the Heaviside step functions with smoothed versions ). The unstable spiral waves of 3V-FK (described by relative equilibria) were found to possess two near-Goldstone modes corresponding to spatial translations with and one Goldstone mode corresponding to temporal translation (or spatial rotation) with . This suggests that translational symmetry is weakly broken due to the presence of boundaries (the calculations were performed on a circular domain of radius equal to just about half the wavelength ). In addition, the spectrum included a pair of complex conjugate modes corresponding to meandering instability and a number of unstable modes corresponding to alternans, all laying on the negative real axis (i.e., corresponding to period-doubling, rather than Hopf, bifurcation). We can therefore expect that the alternans instability may lead to substantially different dynamics in the two models.
The origin of symmetry breaking in the Karma model can be understood by considering the temporal evolution of the group tangents , , and , which become Goldstone modes in the presence of global continuous symmetries. For the Goldstone modes we should have (up to the level of precision defined by ). We checked this by evolving the group tangents using the linearization of (1) about the periodic solution shown in Fig. 3a. As Fig. 4 illustrates, the temporal tangent is a Goldstone mode: while the spatial tangents are not. Both and achieve their maxima on the boundary of (Figs. 4d and 4e), which confirms the role of boundaries in breaking the global translational symmetry. However, there is another mechanism that also breaks translational symmetry. By applying windowing to suppress the boundary effects, we discover that and behave like Goldstone modes everywhere except near the core of the spiral wave (Figs. 4g and 4h), suggesting that translational symmetry is also broken locally. We investigate this local mechanism next.
iv.1 Effect of discretization
As we mentioned previously, spatial discretization can be a local source of symmetry breaking. For the symmetry to get broken the solution should possess a length scale comparable to the scale of spatial discreteness, . The spiral wave solution is characterized by several scales: the wavelength , the length scale of the fast (voltage) variable, which defines the width of the sharp leading front of the excitation wave, and the length scale over which the dynamics of the slow variable switches from excitable to refractory. The latter length scale is controlled by the term in (2) and, to leading order in , is given by . Setting gives . Hence, for () the solution is spatially ill-resolved, and the continuous translational symmetry is broken near the tip, pinning the spiral. On the other hand, for () the solution is well-resolved and continuous translational symmetry should be preserved, so the spiral can drift.
Indeed, this is exactly what we find for . The spiral wave solution at this low value of is similar to that shown in Fig. 3a, but corresponds to a generalized relative periodic orbit with period and wavelength . The displacement of the wave over one period is small, but non-zero (), and its direction depends on both the position and the phase of the spiral wave. The corresponding spectrum (Fig. 5a) contains two complex conjugate unstable eigenvalues which correspond to the alternans instability (the corresponding modes are similar to those shown in Figs. 3e and 3f) and three eigenvalues with . The corresponding Goldstone modes, predictably, coincide with the three group tangents (Fig. 5d), (Fig. 5c), and (Fig. 5b). This indicates that although global Euclidean symmetry is broken on the finite domain, local Euclidean symmetry (i.e., symmetry with respect to small translations or rotations) remains essentially exact provided that (i) the domain is sufficiently large and (ii) the length scale of heterogeneities is sufficiently small.
The transition between the small- regime where local Euclidean symmetry is preserved and the large- regime where it is broken is continuous. The -dependence of the leading eigenvalues is shown in Fig. 6a. As the value of is increased, two of the three unit eigenvalues split off around and separate along the real axis. For this particular solution, Goldstone mode becomes unstable (and frustrated), while the Goldstone mode becomes stable (and frustrated) at large . The exact symmetry of the problem with respect to rotations by means that there is a rotated copy of the solution we found for which the stability of these two modes is interchanged. In either case, however, we find that the and directions play a special role: discretization breaks the rotational symmetry despite the use of the Laplacian stencil that aims to preserve it.
It is also worth pointing out that for , several other modes become unstable. Unlike the modes shown in Figs. 3 and 5 which are isolated (belong to the discrete spectrum), these new modes are not isolated (i.e, belong to the continuous spectrum). Hence, we should expect a continuum of unstable alternans-like modes to appear for single-spiral solutions in unbounded domains for intermediate values of .
The variation of , which defines the net spatial displacement of the wave over one period of rotation, with is shown in Fig. 6b. Consistent with our dimensional analysis, we find that vanishes for , but does not vanish for . Even though the spirals drift for small , when their tip is far from domain boundaries, the displacement over one period is too small to be resolved in DNS and can only be computed using extremely precise calculations based on Newton-Krylov method. However, we will see below that the drift can become quite significant for spirals whose tip is close to a boundary.
The period of the spiral was found to be a monotonically decreasing function of , as Fig. 6c illustrates. The deviation from the limiting value , which corresponds to , was found to be well-approximated by a power law, . The overall variation, however, was not large, with the period being just 12% larger at than at .
To further verify our dimensional analysis we also computed the single-spiral solution on progressively finer grids with fixed stiffness parameter and physical domain size and determined that the solution recovers the Goldstone modes associated with translational symmetry on sufficiently fine grids. The unit eigenvalues are recovered for the translational modes with precision when , whereas the dimensional analysis predicts a similar critical length scale for this value of .
iv.2 Boundary effects
Global translational symmetry implies that there are infinitely many copies of the solution on an unbounded domain that differ in their position but are otherwise equivalent. On bounded domains characterized by local translational symmetry we can also find multiple solutions related to each other by a translation, but they are not, strictly speaking, equivalent. To quantify this relation more precisely we performed a continuation of the domain-centered solutions discussed previously by gradually shifting the tip towards one of the boundaries on a reasonably large domain. The choice of the boundary is arbitrary due to the 4-fold rotational symmetry of the problem.
For the present purposes it is convenient to place the origin of coordinates at the lower left corner of the domain. If the tip of the spiral (identified as the point at which ) has coordinates , then and define the distance of the tip of the spiral, respectively, from the left and bottom boundary. The continuation sequence involves shifting the converged spiral wave solution toward the left or right boundary (decreasing or increasing ) to generate an initial condition, which is refined into a new spiral wave solution using the Newton-Krylov solver. This cycle repeats until the Newton-Krylov solver either fails to converge or converges to a previously found solution.
In the large- regime local continuous symmetry is broken, so we only find solutions shifted by an integer multiple of the grid spacing . Regardless of the direction (towards the left or right wall) the spiral wave solution with shown in Fig. 3a could only be continued for , i.e., no closer than to either wall (cf. Fig. 7c). We will denote this solution branch . When was decreased below about , Newton-Krylov solver converged to a nearby spiral wave with . This new spiral wave solution can be continued until . We will denote this solution branch . If continued in the opposite direction, can be extended symmetrically to .
An example of a spiral wave corresponding to the branch with the tip near along with its spectrum is shown in Fig. 8. Its shape is essentially indistinguishable from that of (cf. Fig. 3a). Its spectrum is also similar to that of (cf. Fig. 3b), except for the eigenvalues with the positive real part: unlike which has a pair of eigenvalues, one unstable and one stable, on the real axis, has two stable complex conjugate eigenvalues . The real and imaginary part of the corresponding modes are shown in Figs. 8e and 8f and can be identified as a linear combination of frustrated translational Goldstone modes (their amplitude varies with distance from the tip). The real and imaginary part of the complex conjugate pair of unstable modes are shown in Figs. 8c and 8d and correspond to the alternans instability.
Counting , , and , there are at least three distinct pinned spiral wave solutions in the large- limit. Their tips are located, modulo , approximately at for , for , and for . Even though these three solutions are distinct, their shapes, temporal periods, and the unstable eigenvalues and eigenmodes corresponding to alternans are virtually indistinguishable, provided they are centered at roughly the same position. Since they can be shifted in either coordinate direction by an integer multiple of , there are “copies” of each of these three solutions.
Most of these “copies” are nearly identical. Consider, for instance the solution . Its leading eigenvalues (cf. Fig. 7a) are effectively independent of the position of the spiral wave and only begin to vary noticeably for . Similarly, the period of this solution (cf. Fig. 7b) is essentially independent of the position of the tip over almost the entire range of . The deviation of the period from the reference value at decreases exponentially fast with : , where is the characteristic length scale that describes the interaction of this spiral with the boundary.
In the small- regime the continuous translational symmetries are recovered, so there is a continuum of drifting spiral wave solutions parametrized by the coordinates of the tip. Unlike the large- limit, there is only one type of solution. Its leading eigenvalues , the spatial shift over one period, and the period as a function of (with ) are shown in Fig. 9. Just like in the large- limit, we find all the basic properties of the spiral wave solution to be effectively independent of the position of the tip, provided . The differences between distinct spiral wave solutions are exponentially small. For instance, Figs. 9b and 9c show that and , where is the period of the centered solution and now . The largest differences ( and ) correspond to the distance of the closest approach .
Our findings can be used to make the concept of local Euclidean symmetry more precise. As long as the tip of the spiral wave is not too close to the boundary of the domain, that solution can be shifted (discretely for large or continuously for small ) without changing any of its properties up to some level of resolution . This level can be made arbitrarily small by increasing the separation between the tip of the spiral waves and the boundary and is only limited (from below) by the numerical tolerance . Hence, for all practical purposes, distinct spiral wave solutions with their tips sufficiently separated from the boundary are completely equivalent, although their spatial shape is affected by the boundaries.
The symmetry breaking in the proximity of the boundary is reflected in the eigenvalues associated with Goldstone modes in the small- limit. As Fig. 9a illustrates, the eigenvalue associated with the -translation deviates from unity for . The vertical boundary does not break the -translation symmetry for spiral waves (and hence does not cause a significant deviation of the corresponding eigenvalue from unity), unless they drift in the direction. This symmetry is eventually broken for when displacement becomes significant. Comparison with Fig. 9a shows that for the corresponding eigenvalue ).
It is worth noting that for , and for , so that the interaction is repulsive (the drift is away from the boundary) at small distances and attractive (the drift is towards the boundary) at slightly larger distances. In a related work Langham and Barkley Langham and Barkley (2013); Langham, Biktasheva, and Barkley (2014) investigated the interaction with the boundaries for resonantly driven (and hence drifting) stable spiral waves in the Barkley model Barkley (1991) of excitable media. They also found that the interaction is repulsive at close range and showed that the interaction length scale is determined by the spatial extent of the response functions (adjoints of the Goldstone modes) which are localized to the core region Biktasheva, Holden, and Biktashev (2006); Biktasheva et al. (2009); Biktasheva, Biktashev, and Foulkes (2010). We have confirmed by direct calculation that both the spatial and temporal response functions in the Karma model (not shown) decay as , where is the distance to the tip.
Finally, the period of rotation shows opposite trends for pinned (large ) and drifting (small ) spirals. As a pinned spiral approaches a boundary it rotates more slowly ( increases), whereas a drifting spiral rotates more quickly ( decreases). A quantitative explanation of these different behaviors relies on the details of the spatial structure of the temporal response function, and is beyond the scope of the current work.
iv.3 Domain size effects
We already have some qualitative intuition about the effect the size of the domain has on the structure and properties of spiral wave solutions. As the size of the domain increases, the solution should approach an unbounded spiral and all its properties should become size-independent. In particular, global Euclidean symmetry (whether continuous for small or discrete for large ) should be restored. On the other hand, as the size of the domain is reduced, the structure of the solution and its properties are expected to change significantly, and, on sufficiently small domains the spiral wave solution might not even exist.
To quantify the differences between single-spiral waves on domains of different size, we computed a sequence of unstable solutions on domains with fixed grid spacing and varying physical size . The sequence was initiated from a large spiral wave with an approximately centered core, which was then truncated on all sides by trimming small regions in a symmetric fashion, which generated the initial condition for the next solution. The process continued until the Newton-Krylov solver failed to find a nontrivial solution.
Local symmetry suggests that sufficiently far from the boundaries the solution should not depend on the boundary condition. Hence, solutions computed on domains of different size, e.g., , should be virtually indistinguishable away from the boundaries of . This is exactly what we find by comparing two solutions and computed on domains of size and , respectively. The difference inside is found to be concentrated in a narrow boundary layer of width (cf. Fig. 10b), where for , as we determined previously. Furthermore, Fig. 10a shows that the difference becomes exponentially small away from the boundaries, , where denotes the distance from the boundary of the smaller domain. Inside the boundary layer the solution on the smaller domain adjusts to the no-flux condition and can have a large curvature comparable to that at the spiral wave core. Hence, it should not be surprising to see the same length scale describe both the width of the core and the width of the boundary layer.
For , the spiral waves are pinned, and we find a minimal domain size, , approximately cm in dimensional units, below which spiral wave solutions cannot be found. This domain size corresponds to the distance between the tip of the spiral and the boundary equal to . Since determines both the radius of the spiral wave core and the width of the boundary layer, corresponds to the collision of the spiral core with the boundary layer. If this criterion applies more generally, we should expect to find spiral waves solutions only on domains of size .
Pinned spiral waves remain unstable in the entire range of system sizes. The dominant eigenvalues of the solution are shown in Fig. 11a. As increases, the eigenvalues from the discrete part of the spectrum approach a constant value, but new eigenvalues also appear which correspond to the continuous part of the spectrum. For , there are between three and five unstable modes. As , the dominant eigenvalues quickly grow in magnitude, approaching (these are outside of the range of Fig. 11a).
Interestingly, the character of the instability changes for . Consider, for instance, the solution at shown in Fig. 12a. The magnitude of its unstable eigenvalues (cf. Fig. 12b) is comparable to that in much larger systems, but the corresponding eigenmodes (Fig. 12c and 12d) correspond to the meandering instability, rather than alternans. Initial conditions close to this unstable solution produce spiral waves which persist for up to rotations without breaking up, which is consistent with numerical results of Karma Karma (1994) for domains smaller than cm (which corresponds to in nondimensional units). As. Fig. 13 illustrates, the amplitude of meandering slowly grows, until the tip runs into a boundary and the wave eventually collapses. The minimal domain size that supports spiral wave breakup via alternans defines another important dynamical length scale.
The period of the solution approaches the asymptotic value exponentially fast as increases (cf. Fig. 11b). For we find , which is consistent with the scaling we found previously. Indeed, for a domain-centered spiral wave , or , which means that scaling with and follows from the same general scaling law. For , the period can be considered essentially independent of .
Most of the results for small are qualitatively similar, so we will only focus on what is new or different. This regime is characterized by continuous local symmetries, with drifting spiral waves possessing three Goldstone modes. The magnitude of the spatial displacement of the spiral wave over one period decreases exponentially with , (cf. Fig. 14b). As the domain is truncated, the eigenvalues associated with Goldstone modes begin to deviate from unity (cf. Fig. 14a). Eigenvalues corresponding to translational modes deviate first, at . This agrees well with the distance to the boundaries at which translational symmetry is broken, as we found in the previous section. Somewhat unexpectedly, for even smaller domains (), the eigenvalue associated with temporal translation also deviates from unity. The critical domain size matches the distance of the tip to the boundary () beyond which the spiral solution loses translational symmetry, as noted in the previous section. Hence deviation of the “temporal” eigenvalue from unity can be understood as a result of the non-perturbative nature of finite spatial shifts on small domains and the subsequent loss of temporal periodicity. For the displacement of the spiral is significant enough to affect its temporal dynamics. The solutions in this limit do not describe relative periodic orbits. For instance, both the “period” and the spatial shift of and becomes noticeably different.
The decrease in the temporal period at small separations between the tip of the spiral and the boundaries observed in Figs. 9c, 11b, and 14c is due to the negative curvature of the wave front near the boundary due to the no-flux boundary condition. It is well-known Keener (1986) that the speed of propagation of an excitation wave is related to the curvature of its front by a linear relationship, with convex () wave fronts traveling more slowly than flat () ones, and flat wave fronts traveling more slowly than concave () ones Zykov and Morozova (1979); Fast and Kléber (1997). The magnitude of this effect on the rotation period is controlled by the spatial structure of the temporal response function, which decays exponentially fast with the distance from the tip of the spiral, with the decay rate equal to the length scale .
The exponential dependence of the shift madnitude and the period of the drifting spiral waves found for the modified Karma model are in contradiction with the analytical results obtained by Aranson et. al. Aranson, Kessler, and Mitkov (1995) for a similar model in the limit, which predict a super-exponential dependence of the drift speed and rotational frequency of spiral waves on the domain size . Davydov and Zykov Davydov and Zykov (1993) predict the frequency of spirals in a generic reaction-diffusion model with to vary as the inverse of the domain size (for circular domains of radius comparable to ). Hartmann et. al. Hartmann et al. (1996) claim that their simulations confirm this prediction, but a quick inspection of their numerical results, as well as those of Davydov and Zykov, shows that their data is in excellent agreement with the exponential dependence.
V Summary and discussion
We performed the first systematic investigation of unstable spiral wave solutions of a simple spatially discretized model of cardiac tissue with physiologically and dynamically relevant no-flux boundary conditions. We also characterized how the basic properties of these solutions, such as their temporal period, spatial drift, and stability, depend on the size of the domain, the proximity of the spiral core to the nearest domain boundary, and the microscopic heterogeneity associated with spatial discretization. We found that, although both the boundary conditions and the discretization break the global Euclidean symmetry of the model, unstable spiral wave solutions tend to respect a local – continuous or discrete – version of Euclidean symmetry.
Existing tools, such as Newton-Krylov solvers, developed for computing unstable solutions in the presence of global continuous symmetries (e.g., in unbounded domains or domains with periodic boundary conditions) break down on bounded domains with generic boundary conditions that are not consistent with the global symmetries. However, for reaction-diffusion systems in general, and monodomain models of cardiac tissue in particular, which lack long-range correlations, the effect of boundaries is localized, which enables computation of solutions satisfying local continuous symmetries using a generalization of Newton-Krylov solvers. The generalization involves weighting, or windowing, of the residual to suppress the symmetry-breaking effect of the boundaries. The generalized solvers permit the computation of unstable solutions describing, e.g., pinned and drifting spiral waves in a model of cardiac tissue in the regime characterized by spontaneous breakup of spiral waves leading to fibrillation-like dynamics.
The inherent spatial heterogeneity associated with the cellular structure of cardiac tissue, and the associated discreteness of the numerical model, was found to have an interesting effect which is usually overlooked in the studies of stable solutions. Cardiac tissue models typically involve discontinuous functions which describe switching of the state of various ionic pumps and channels. Although they simplify the formulation, these discontinuities are unphysical and lead to the emergence of short time scales and corresponding short length scales that can impact the properties of unstable solutions. When these short length scales are small compared to the size of cardiomyocytes, the spatial heterogeneity associated with cellular structure breaks the continuous translational symmetry, leading to pinning of unstable spiral waves. In the particular cardiac model considered here, at least three different branches of such pinned spiral waves exist, parametrized by the position of their tip relative to the computational grid and characterized by distinct stability properties. Notably, the “alternans-stable” meandering spiral waves found in the same parameter regime do not display pinning. However, even there discreteness manifests itself in the jaggedness of the tip trajectory.
Even though microscopic spatial heterogeneity breaks the continuous translational symmetry, discrete translation symmetry survives on sufficiently large domains. We find that each unstable solution can be shifted by a discrete number of grid points (or cells) along or transversely to the direction of the fibers without changing, to numerical precision, either its stability properties or its period. Hence, for all practical purposes these solutions can be considered equivalent under discrete translations. This translation symmetry is local in the sense that the properties of different solutions from the same branch are only invariant (to numerical precision) for finite translations such that the tip of the spiral wave remains outside of the boundary layer of width , but they begin to vary as the tip approaches any of the boundaries. For each branch, the translation symmetry is reflected in the stability spectrum in the form of slightly frustrated analogues of translational Goldstone modes characterized by a pair of real or complex conjugate Floquet multipliers with positive real part.
Continuous translational symmetry can be restored by replacing the discontinuous switching functions in the ionic model with smooth ones. In the absence of small intrinsic length scales, spiral wave solutions become well-resolved even on a discrete grid. As a result, three discrete branches of pinned spiral wave solutions are merged into a single branch which comprises a continuum of symmetry-related drifting spiral waves parametrized by the position of their tip. Here too, the continuous translational symmetry is manifested in the stability spectra of the solutions, which possess three Goldstone modes corresponding to the three continuous translation symmetries, one with respect to time and two with respect to spatial position. Again, the spatial symmetry is local: the properties of spiral wave solutions remain invariant (to numerical precision) with respect to finite translations, provided the tip of the spiral wave remains outside of the boundary layer of width .
Global Euclidean symmetry can, of course, be gradually restored by increasing the size of the computational domain. However, even on finite domains one finds that the solutions approach their asymptotic shape in the interior of the domain exponentially fast as increases. The difference between solutions with different is again found to be confined to the boundary layer of width . Similarly, the properties of solutions approach the asymptotic values exponentially fast as increases, with the same length scale . This length scale () was found to control the effect of the boundaries rather generally. Effectively, controls the strength of interaction between a spiral wave and a (locally) straight boundary (and by extension, interaction between two identical counter-rotating spiral waves). With few exceptions (notably, the eigenvalues from the continuous part of the stability spectrum), the deviation from asymptotic values for all properties of spiral wave solutions (e.g., their shapes, temporal periods, spatial displacement over one rotation) were found to scale exponentially, i.e., as , with the distance to the boundary.
In conclusion, let us comment on the implications of our results for the problem of fibrillation. Unstable spiral wave solutions can only be found for sufficiently large domains, , with the smallest domain size corresponding to the strongly-interacting regime in the Karma model. On domains smaller than this size, spiral waves do not persist for a complete rotation. This length scale determines, in both limits of the stiffness parameter considered here, the minimal spacing between spiral cores in the state of fibrillation Byrne, Marcotte, and Grigoriev (2015). As varies between and , the period of the spiral wave varies as much as 20%, with the smaller spirals rotating faster than the larger ones. Hence, even if we were to ignore the instabilities on time scales of order a few periods, multi-spiral states with spirals of different size should exhibit dynamics that are very complicated, i.e., at least quasiperiodic. Quasiperiodicity has been suggested as a possible dynamical mechanism for transition to fibrillation Garfinkel et al. (1997).
We have established that, for meander is the leading instability mechanism, so spirals of size will tend to merge with neighboring spirals of opposite chirality in the same size range. This mechanism reduces the number of spirals and increases the sizes of remaining spirals. On the other hand, spirals of size are unstable towards alternans and will break up. This mechanism increases the number of spirals by splitting large spirals into smaller ones (ranging in size down to ). The interaction of these two instability mechanisms alone can lead to a dynamic self-sustaining process, which would maintain the state of spiral chaos featuring multiple interacting spirals ranging in size from to .
Acknowledgements.
The authors would like to thank Vadim Biktashev and Dwight Barkley for helpful discussions and comments. This material is based upon work supported by the National Science Foundation under Grant No. CMMI-1028133. The Tesla K20 GPUs used for this research were donated by the “NVIDIA Corporation” through the academic hardware donation program.Vi References
References
- Cherry and Fenton (2008) E. M. Cherry and F. H. Fenton, “Visualization of spiral and scroll waves in simulated and experimental cardiac tissue,” New J. Phys. 10, 125016 (2008).
- Gray, Pertsov, and Jalife (1998) R. A. Gray, A. M. Pertsov, and J. Jalife, “Spatial and temporal organization during cardiac fibrillation,” Nature 392, 75–78 (1998).
- Clayton et al. (2011) R. H. Clayton, O. Bernus, E. M. Cherry, H. Dierckx, F. H. Fenton, L. Mirabella, A. V. Panfilov, F. B. Sachse, G. Seemann, and H. Zhang, “Models of cardiac tissue electrophysiology: Progress, challenges and open questions,” Prog. Biophys. Mol. Biol. 104, 22–48 (2011).
- Courtemanche (1996) M. Courtemanche, “Complex spiral wave dynamics in a spatially distributed ionic model of cardiac electrical activity,” Chaos 6, 579 (1996).
- Barkley (1991) D. Barkley, “A model for fast computer simulation of waves in excitable media,” Physica D 49, 61–70 (1991).
- Karma (1994) A. Karma, “Electrical alternans and spiral wave breakup in cardiac tissue,” Chaos 4, 461–472 (1994).
- Bueno-Orovio, Cherry, and Fenton (2008) A. Bueno-Orovio, E. M. Cherry, and F. H. Fenton, ‘‘Minimal model for human ventricular action potentials in tissue,” J. Theor. Biol. 253, 544–560 (2008).
- Simitev and Biktashev (2006) R. D. Simitev and V. N. Biktashev, “Conditions for propagation and block of excitation in an asymptotic model of atrial tissue,” Biophysical journal 90, 2258–2269 (2006).
- Noble (1962) D. Noble, “A modification of the Hodgkin-Huxley equations applicable to Purkinje fibre action and pace-maker potentials,” J. Physiol. 160, 317–352 (1962).
- Cherry et al. (2007) E. M. Cherry, J. R. Ehrlich, S. Nattel, and F. H. Fenton, “Pulmonary vein reentry-Properties and size matter: Insights from a computational analysis,” Heart Rhythm 4, 1553–1562 (2007).
- Muller, Plesser, and Hess (1985) S. C. Muller, T. Plesser, and B. Hess, “The structure of the core of the spiral wave in the Belousov-Zhabotinskii reaction,” Science 230, 661–663 (1985).
- Keener and Tyson (1986) J. P. Keener and J. J. Tyson, “Spiral waves in the Belousov-Zhabotinskii reaction,” Physica D 21, 307–324 (1986).
- Siegert and Weijer (1991) F. Siegert and C. J. Weijer, “Analysis of optical density wave propagation and cell movement in the cellular slime mold Dictyostelium discoideum,” Physica D 49, 224–232 (1991).
- Vasiev, Hogeweg, and Panfilov (1994) B. N. Vasiev, P. Hogeweg, and A. V. Panfilov, “Simulation of Dictyostelium discoideum aggregation via reaction-diffusion model,” Phys. Rev. Lett. 73, 3173–3176 (1994).
- Murray, Stanley, and Brown (1986) J. D. Murray, E. A. Stanley, and D. L. Brown, “On the spatial spread of rabies among foxes,” Proc. Roy. Soc. London. B 229, 111–150 (1986).
- Ouyang and Swinney (1991) Q. Ouyang and H. L. Swinney, “Transition to chemical turbulence,” Chaos 1 (1991), 10.1063/1.165851.
- Hagberg and Meron (1994) A. Hagberg and E. Meron, “From labyrinthine patterns to spiral turbulence,” Phys. Rev. Lett. 72, 2494–2497 (1994).
- Panfilov and Hogeweg (1993) A. V. Panfilov and P. Hogeweg, “Spiral breakup in a modified FitzHugh-Nagumo model,” Phys. Lett. A 176, 295–299 (1993).
- Morris et al. (1993) S. W. Morris, E. Bodenschatz, D. S. Cannell, and G. Ahlers, “Spiral defect chaos in large aspect ratio Rayleigh-Benard convection,” Phys. Rev. Lett. 71, 2026–2029 (1993).
- Barkley, Kness, and Tuckerman (1990) D. Barkley, M. Kness, and L. S. Tuckerman, “Spiral wave dynamics in a simple model of excitable media: Transition from simple to compound rotation,” Phys. Rev. A 42, 2489–2492 (1990).
- Barkley (1992) D. Barkley, “Linear stability analysis of rotating spiral waves in excitable media,” Phys. Rev. Lett. 68, 2090–2093 (1992).
- Goldhirsch, Orszag, and Maulik (1987) I. Goldhirsch, S. A. Orszag, and B. K. Maulik, “An efficient method for computing leading eigenvalues and eigenvectors of large asymmetric matrices,” J. Sci. Comp. 2, 33–58 (1987).
- Allexandre and Otani (2004) D. Allexandre and N. F. Otani, “Preventing alternans-induced spiral wave breakup in cardiac tissue: An ion-channel-based approach,” Phys. Rev. E 70, 061903 (2004).
- Barkley (1994) D. Barkley, “Euclidean symmetry and the dynamics of rotating spiral waves,” Phys. Rev. Lett. 72, 164–167 (1994).
- Barkley and Kevrekidis (1994) D. Barkley and I. G. Kevrekidis, “A dynamical systems approach to spiral wave dynamics,” Chaos 4, 453–460 (1994).
- Fiedler et al. (1996) B. Fiedler, B. Sandstede, A. Scheel, and C. Wulff, “Bifurcation from relative equilibria of noncompact group actions: Skew products, meanders, and drifts,” Doc. Math. 141, 479–505 (1996).
- Fiedler and Turaev (1998) B. Fiedler and D. Turaev, “Normal forms, resonances, and meandering tip motions near relative equilibria of Euclidean group actions,” Arch. Rational Mech. Anal. 145, 129–159 (1998).
- Sandstede, Scheel, and Wulff (1997) B. Sandstede, A. Scheel, and C. Wulff, “Dynamics of spiral waves on unbounded domains using center-manifold reductions,” J. Diff. Eqn. 141, 122–149 (1997).
- Sandstede, Scheel, and Wulff (1999) B. Sandstede, A. Scheel, and C. Wulff, “Dynamical behavior of patterns with Euclidean symmetry,” in Pattern Formation in Continuous and Coupled Systems (Springer, New York, 1999) pp. 249–264.
- Chossat and Lauterbach (2000) P. Chossat and R. Lauterbach, Methods in Equivariant Bifurcations and Dynamical Systems (World Scientific, Singapore, 2000).
- Beyn and Thümmler (2004) W.-J. Beyn and V. Thümmler, “Freezing solutions of equivariant evolution equations,” SIAM J. Appl. Dyn. Syst. 3, 85–116 (2004).
- Beyn and Lorenz (2008) W.-J. Beyn and J. Lorenz, “Nonlinear stability of rotating patterns,” Dynamics of PDE 5, 349–400 (2008).
- Hermann and Gottwald (2010) S. Hermann and G. A. Gottwald, “The large core limit of spiral waves in excitable media: A numerical approach,” SIAM J. Appl. Dyn. Sys. 9, 536–567 (2010), arXiv:1003.5830.
- Foulkes and Biktashev (2010) A. J. Foulkes and V. N. Biktashev, ‘‘Riding a spiral wave: Numerical simulation of spiral waves in a comoving frame of reference,” Phys. Rev. E 81, 046702 (2010), arXiv:1001.4454.
- Fenton, Cherry, and Glass (2008) F. H. Fenton, E. M. Cherry, and L. Glass, “Cardiac arrhythmia,” Scholarpedia 3, 1665 (2008).
- Bär and Eiswirth (1993) M. Bär and M. Eiswirth, “Turbulence due to spiral breakup in a continuous excitable medium,” Phys. Rev. E 48, R1635–R1637 (1993).
- Fenton et al. (2002) F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J. Evans, “Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity,” Chaos 12, 852–892 (2002).
- Bernus, Verschelde, and Panfilov (2003) O. Bernus, H. Verschelde, and A. V. Panfilov, “Spiral wave stability in cardiac tissue with biphasic restitution,” Physical Review E 68, 021917 (2003).
- Biktashev (2007) V. N. Biktashev, “Drift of spiral waves,” Scholarpedia (2007).
- Fenton and Cherry (2008) F. H. Fenton and E. M. Cherry, “Models of cardiac cell,” Scholarpedia 3, 1868 (2008).
- Sepulveda, Roth, and Jr. (1989) N. G. Sepulveda, B. J. Roth, and J. P. W. Jr., “Current injection into a two-dimensional anisotropic bidomain,” Biophys. J. 55, 987–999 (1989).
- Roth (2008) B. J. Roth, “Bidomain model,” Scholarpedia 3, 6221 (2008).
- Karma (1993) A. Karma, “Spiral breakup in model equations of action-potential propagation in cardiac tissue,” Phys. Rev. Lett. 71, 1103–1106 (1993).
- Bevans et al. (1998) C. G. Bevans, M. Kordel, S. K. Rhee, and A. L. Harris, “Isoform composition of connexin channels determines selectivity among second messengers and uncharged molecules,” J. Biol. Chem. 273, 2808–2816 (1998).
- Garcia-Dorado, Rodriguez-Sinovas, and Ruiz-Meana (2004) D. Garcia-Dorado, A. Rodriguez-Sinovas, and M. Ruiz-Meana, “Gap junction-mediated spread of cell injury and death during myocardial ischemia-reperfusion,” Cardiovasc. Res. 61, 386–401 (2004).
- Guevara et al. (1984) M. R. Guevara, G. Ward, A. Shrier, and L. Glass, “Electrical alternans and period-doubling bifurcations,” Comput. Cardiol. , 167 (1984).
- Nolasco and Dahlen (1968) J. B. Nolasco and R. W. Dahlen, “A graphic method for the study of alternation in cardiac action potentials,” J. Appl. Physiol. 25, 191–196 (1968).
- Byrne, Marcotte, and Grigoriev (2015) G. Byrne, C. D. Marcotte, and R. O. Grigoriev, “Exact coherent structures and chaotic dynamics in a model of cardiac tissue,” Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 033108 (2015).
- Biktashev, Holden, and Nikolaev (1996) V. N. Biktashev, A. V. Holden, and E. V. Nikolaev, “Spiral wave meander and symmetry of the plane,” Int. J. Bifur. Chaos 6, 2433–2440 (1996).
- Biktasheva et al. (2009) I. V. Biktasheva, D. Barkley, V. N. Biktashev, G. V. Bordyugov, and A. J. Foulkes, “Computation of the response functions of spiral waves in active media,” Phys. Rev. E 79, 056702 (2009).
- Biktasheva, Biktashev, and Foulkes (2010) I. V. Biktasheva, V. N. Biktashev, and A. J. Foulkes, “Computation of the drift velocity of spiral waves using response functions,” Phys. Rev. E 81, 066202 (2010).
- Lindeberg (1990) T. Lindeberg, “Scale-space for discrete signals,” IEEE Transactions on Pattern Analysis and Machine Intelligence 12, 234–254 (1990).
- Viswanath (2007) D. Viswanath, “Recurrent motions within plane Couette turbulence,” J. Fluid Mech. 580, 339–358 (2007), arXiv:physics/0604062.
- Knoll and Keyes (2004) D. A. Knoll and D. E. Keyes, “Jacobian-free Newton-Krylov methods: A survey of approaches and applications,” J. Comp. Phys. 193, 357 (2004).
- Bueno-Orovio, Perez-Garcia, and Fenton (2006) A. Bueno-Orovio, V. M. Perez-Garcia, and F. H. Fenton, “Spectral methods for partial differential equations in irregular domains: The spectral smoothed boundary method,” SIAM J. Sci. Comp. 28, 886–900 (2006).
- Bjorstad and Gropp (2004) P. Bjorstad and W. Gropp, Domain decomposition: parallel multilevel methods for elliptic partial differential equations (Cambridge Univ. Press, 2004).
- Xu (1992) J. Xu, “Iterative methods by space decomposition and subspace correction,” SIAM review 34, 581–613 (1992).
- Teramura and Toh (2014) T. Teramura and S. Toh, “Damping filter method for obtaining spatially localized solutions,” Physical Review E 89, 052910 (2014).
- Trefethen and Bau (1997) L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM, 1997).
- Fenton and Karma (1998) F. Fenton and A. Karma, “Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation,” Chaos 8, 20–47 (1998).
- Langham and Barkley (2013) J. Langham and D. Barkley, “Non-specular reflections in a macroscopic system with wave-particle duality: Spiral waves in bounded media,” Chaos 23, 013134 (2013), arXiv:1304.0591.
- Langham, Biktasheva, and Barkley (2014) J. Langham, I. Biktasheva, and D. Barkley, “Asymptotic theory for spiral wave reflections,” (2014), arXiv:11401.7626.
- Biktasheva, Holden, and Biktashev (2006) I. V. Biktasheva, A. V. Holden, and V. N. Biktashev, “Localization of response functions of spiral waves in the FitzHugh–Nagumo system,” Int. J. Bifur. Chaos 16, 1547–1555 (2006).
- Keener (1986) J. P. Keener, “A geometrical theory for spiral waves in excitable media,” SIAM J. Appl. Math. 46, 1039–1056 (1986).
- Zykov and Morozova (1979) V. S. Zykov and O. L. Morozova, “Speed of spread of excitation in a two-dimensional excitable medium,” Biophysics 24, 739–744 (1979).
- Fast and Kléber (1997) V. G. Fast and A. G. Kléber, “Role of wavefront curvature in propagation of cardiac impulse,” Cardiovascular research 33, 258–271 (1997).
- Aranson, Kessler, and Mitkov (1995) I. Aranson, D. Kessler, and I. Mitkov, “Drift of spiral waves in excitable media,” Physica D: Nonlinear Phenomena 85, 142–155 (1995).
- Davydov and Zykov (1993) V. A. Davydov and V. S. Zykov, “Spiral autowaves in a round excitable medium,” J. Exp. Theor. Phys. 76, 414–419 (1993).
- Hartmann et al. (1996) N. Hartmann, M. Bär, I. Kevrekidis, K. Krischer, and R. Imbihl, “Rotating chemical waves in small circular domains,” Physical review letters 76, 1384 (1996).
- Garfinkel et al. (1997) A. Garfinkel, P. S. Chen, D. O. Walter, H. S. Karagueuzian, B. Kogan, S. J. Evans, M. Karpoukhin, C. Hwang, T. Uchida, M. Gotoh, O. Nwasokwa, P. Sager, and J. N. Weiss, “Quasiperiodicity and chaos in cardiac fibrillation,” J. Clin. Invest. 99, 305–314 (1997).