Unstable spiral waves and local Euclidean symmetry in a model of cardiac tissue

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 coherences of the contractile dynamics are destroyed. Ventricular fibrillation^1,2 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.^3–10 Besides cardiac tissue, reaction-diffusion systems are also used to model diverse phenomena such as chemical reactions,^11,12 bacterial chemotaxis,^13,14 and disease propagation.¹⁵ 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,¹⁶ spiral chaos,¹⁷ spiral breakup,¹⁸ or spiral defect chaos.¹⁹ 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.²⁰ 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²¹ 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²² 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.²³ 

Barkley²⁴ and Barkley and Kevrekidis²⁵ showed that the tip dynamics for both pinned and meandering spirals can be reproduced using a system of five weakly nonlinear ordinary differential equations (ODEs) derived by assuming the dynamics are equivariant under a Euclidean symmetry group. Fiedler et al.^26,27 showed that, more generally, equivariance with respect to noncompact, finite-dimensional Lie groups (such as the Euclidean group E(n)) 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²⁸ and near relative periodic orbits,²⁹ formalizing and extending the reduced description of Barkley and Kevrekidis^24,25 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.³⁰ 

Beyn and Thümmler³¹ 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³² 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³³ to investigate the dynamics of spiral waves in the large-core limit and by Foulkes and Biktashev³⁴ 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³⁵ and as a result are typically strongly unstable, often encountering refractory regions of tissue and breaking up within a few rotations.^36–38 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 Sec. 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 Sec. 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 Sec. IV. Section V summarizes these results and discusses their implications in the context of fibrillation (or, more generally, spiral chaos).

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³⁹ 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 n interacting chemicals, u = [u₁, u₂, … ,u_n], their evolution can be described by the system of partial differential equations (PDEs)

where D 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, u = [u, v], where u is an (electric) transmembrane potential and v 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.⁴⁰ Physically more realistic bidomain models⁴¹ 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.⁴² 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.⁶ 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.^6,43 We use the following nondimensionalized form of the “reaction” terms:

where there is only one gating variable, so that u = [u, v].

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,²³ with a smooth function $Θs(u)=(1+tanh(su))/2$⁠. This substitution removes the singularities in the partial derivatives ∂f_i/∂u_j, while recovering the original step function in the limit s →∞. An additional term Θ_s(v − 1)(v − 1) was introduced to smooth out an unphysical singularity in the v 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 D₂₂ = νD₁₁ with 0 < ν ≪ 1, 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,^44,45 so a small nonzero value of D₂₂ is justified from the physiological perspective. These modifications have a very subtle effect on the cellular dynamics. In particular, the v-nullcline (Fig. 1(a)) is smoothed into a sharply varying, but C^∞-continuous function (Fig. 1(b)). The position and the stability of the equilibria remain essentially unchanged for s sufficiently large, however, significant deviation in the position of the equilibria is found as s is decreased to O(1).

We kept the parameter values used in the original work,⁶ i.e., u* = 1.5415, M = 4, and ϵ = 0.01. All times and lengths were non-dimensionalized using, respectively, the characteristic time scale τ_u = 2.5 ms and the characteristic length scale ξ = 0.0262 cm (the size of a typical cardiomyocyte). The original diffusion constant γ = 1.1 cm²/s was non-dimensionalized accordingly, yielding D₁₁ = γτ_u/ξ²= 4.0062. The ratio of the two diffusion constants was set to ν = 0.05. Finally, the width s^–1 of the smoothed step function Θ_s(u) was varied in the range 1 ≤ s ≤ 32.

The dynamics of the modified model with s ≳ 4 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 $R=ln(β/(β−1))$⁠. Karma showed⁶ that the pinned spiral solution is stable at intermediate values of R. As R 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 R 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.^46,47 As R 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 R = 1.273 (which corresponds to β = 1.389). For this value of R, spiral chaos persists indefinitely for sufficiently large domains.⁶ 

We used the physiologically and dynamically relevant no-flux boundary conditions, n⋅∇u = 0. 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⁴⁸ of multi-spiral states characterizing spiral chaos showed that n⋅∇u ≈ 0 at the edges of domains (or tiles) supporting individual spirals.

Reaction-diffusion systems (1) with arbitrary local kinetics are equivariant with respect to Euclidean transformations $g∈G=G×E+(1)$⁠, where $G=E(n)$ for an n-dimensional physical space, i.e., $ℝn$⁠. The subgroup E⁺(1) reflects the invariance with respect to time translation, but not time-reversal. In an unbounded two-dimensional domain, the corresponding symmetry group E(2) × E⁺(1) includes continuous translations along the two coordinate directions (x, y), 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 ∂_tu − ω∂_θu = 0, where ω is the angular speed and ∂_θ = (x − x₀)∂_y − (y − y₀)∂_x is the angular derivative about an origin [x₀, y₀]. 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 T, gu(T) = u(0), up to a symmetry transformation $g∈G$⁠. 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 ℓ_c of its core, which is determined by the spatial extent of the response functions (adjoint eigenfunctions with unit Floquet multipliers).^49–51 In particular, if the distance between the tip of the spiral wave and the boundary is larger than O(ℓ_c), 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 Sec. III using square domains as a fairly representative example.

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 x (y) axis is oriented along (transversely to) the fibers.

The modified Karma model (1) and (2) is solved numerically on square domains of side length L by using a finite-difference discretization of the Laplacian operator on a two-dimensional grid of size N × N. Since we chose the size of cardiac cells as the characteristic length scale, ξ, by construction the grid spacing is Δx = Δy = 1 and L = N 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 u at position (x, y) is approximated using a second-order finite-difference stencil

with coefficients a_±1,0 = a_0,±1 = ζ, a_±1,±1 = (1 − ζ)/2, and a_0,0 = −2(1 + ζ). For ζ = 2/3, this expression is the most isotropic formulation of the nine-point finite-difference stencil on a two-dimensional uniform grid,⁵² 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 Δt = 0.004. This time step is sufficiently small to accurately resolve the fast time scales in the evolution of the u field, and found to produce solutions accurate to one part in 10¹⁰ in the 2-norm for integration times of order 10⁴Δt (or 100 ms, in dimensional units) with the expected O(Δt⁴) 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.^21,23 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⁵³ 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.

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 T, but shifted by $h=x̂hx+ŷhy$ relative to the initial condition (generalization to the case of meandering spiral waves involving rotation is straightforward). Using the evolution operator $UT=exp(T∂t)$ and the operator of spatial translations $Th=exp(h·∇)$⁠, this can be written as $T−hUTu=u$⁠. This solution can be specified completely by a real-valued vector $w=[u,hx,hy,T]$⁠, where the state $u=[u1,…,uN2,v1,…,vN2]$ is composed of the values of the fields u and v on the computational grid. More generally, a point that lies on a relative periodic orbit corresponds to a root of the vector function $F:ℝ2N2×ℝ|G|→ℝ2N2$⁠,

which measures the difference between the “initial” state u and the transformed “final” state $g UTu$ (here, the more general group transformation $g∈G$ takes the place of the translation operator $T−h$⁠). Absolute (non-relative) periodic solutions are merely a special case: they satisfy Eq. (4) with $g=1$⁠, (i.e., h = 0). Similarly, relative equilibria may be treated as special cases of periodic orbits.

Taylor expansion of Eq. (4) shows that the leading-order approximation for the correction δu to the system state is given by the solution of

where $A=∂uF|w$⁠. However, Eq. (5) does not determine the $|G|$ corrections specifying the group transformation and the periodicity, (e.g., the δh and δT of our drifting spiral). This is a generic feature of systems with continuous symmetries, since there is a continuum of solutions satisfying Eq. (4). Additional constraints are needed to uniquely specify the solution. Specifically, we will require that the correction to the state δu be transverse to the group manifold. The constraints can be written using the generators of translation along x, y, and t, which are given by the respective partial derivatives, yielding the following well-posed linear problem:

where $F′=[F,0]$⁠, and

where $u′=g UTu$ denotes the final state, and $∇q$ is the row vector of partial derivatives with respect to the independent coordinates q which define the local tangents of the group manifold, e.g., $∇q=(∂x,∂y,∂t)$ for the drifting spiral.

The linear system (6) effectively implements the skew-product representation of the dynamics, where the correction $δw=[δu,δq]$ is split into components transverse to, and along, the group manifold. Many authors⁵⁴ drop the dependence on the final state altogether, instead substituting $∇qu′→∇qu$⁠. We are unaware of any benefit in doing so and have opted not to, as $u′$ is explicitly available: $u′=u+F$⁠. 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.²¹ For example, in the special case of relative equilibria we have ∂_θu = ω⁻¹∂_t u.

Typical discretizations of Eq. (2) involve 10⁴–10⁶ variables, making it very expensive to even compute the elements of the matrix $A′$⁠, let alone solve 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,²² 

where the action of A is approximated using finite-differencing with some $|α|≪1$ and $v̂i$ denotes the normalized projection of v_i onto the orthogonal complement of the linear space spanned by ${v̂1,…,v̂i−1}$⁠.

Projecting Eq. (6) into the Krylov subspace spanned by ${v̂1,…,v̂k}$⁠, we obtain a $(k+|G|)$-dimensional linear system

where $Hk=PkAPk†$ is a real-valued k × k upper-Hessenberg matrix generated explicitly during the Arnoldi iteration, $cq=Pk∇qu, aq=Pk∇qu′, b=PkF, δy=Pkδu$⁠, and $Pk=[v̂1⋯v̂k]†$ is the projection operator. The vector δq = [δh, δT] defines the shift and period corrections in our drifting spiral example. Due to strong contraction associated with the Laplacian term in Eq. (1), accurate approximation of the original system (6) can usually be obtained with k of order a few tens, very small compared with the full dimensionality of the discretized system, yielding the solution $δw≈[Pk†δy,δq]$⁠.

Newton-Krylov iterations w_n+1 = w_n + δw_n converge provided a good initial condition w₀ 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, $∥F(wn+δwn)∥>∥F(wn)+A δwn∥$⁠. A scaling factor 0 < η_n ≤ 1 for the Newton step magnitude was computed by minimizing the 2-norm of the residual $F(wn+ηnδwn)$⁠, yielding a more robust sequence of iterations $wn+1=wn+ηnδwn$⁠. The iteration is considered converged to the solution w when the 2-norm of the residual function is minimized below a predetermined value, $‖F(w)‖<εtol$ (in this study, we used ε_tol = 10⁻¹⁰).

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 h = δh = 0 and dropping the condition of transversality with respect to translations, so that δq → δT. However, Newton-Krylov iterations predictably fail to converge for drifting spirals for which h ≠ 0. Since boundary conditions break translational symmetry, the residual (4) cannot be reduced below $C|h|$ for some positive constant C.

Boundary conditions do not merely break the translational and rotational symmetries. Finite rotations $Rϕ=exp(ϕ∂θ)$ on a bounded domain Ω are not injective: some points are mapped out of the domain and others into (Fig. 2(a)). A similar situation is encountered for finite translations $T−h$ on a bounded domain (Fig. 2(b)). Consider, for example, a meandering spiral wave for which $RϕUTu=u$ on an unbounded domain. On a bounded domain, the residual $F=RϕUTu−u$ will not vanish (we set u ≡ 0 outside of Ω to make the residual well-defined). If one places the origin of rotation near the tip of the spiral wave, the residual F will decompose into two easily identifiable contributions. Inside $Ω∩RϕΩ$ (the octagonal overlap region in Fig. 2(a)), the residual is small and lies near the group manifold,

where δq 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 $(Ω∪RϕΩ)∖(Ω∩RϕΩ)$ (the eight triangular regions in Fig. 2(a)), the residual is large, F = O(1).

Figure 2(b) describes the effect of boundaries on drifting spiral waves described by relative periodic orbits for which $T−hUTu=u$ on an unbounded domain. On a bounded domain, we find a decomposition of the residual $F=T−hUTu−u$ analogous to the case of meandering waves. In the overlap region $Ω∩ThΩ$ (the rectangular region in Fig. 2(b)), 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, F = O(1). 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,

where ρ is a diagonal matrix with elements 0 ≤ ρ_ii ≤ 1. The first 2N² diagonal elements correspond to the weights associated with the dynamical field variables u and v; they correspond to the values of the windowing function

where x_c denotes the center of the domain and d determines the diameter of the (circular) “window” in units of L (we set d = 0.7 and σ = 32, unless specified otherwise). Specifically, ρ_ii = W_d(x_i). 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 Figs. 2(c) and 2(d), 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,⁵⁵ domain decomposition,^56,57 and the damping filter method.⁵⁸ Alternatively, the weighting may be interpreted as an ad hoc Jacobi-type preconditioning method for solving linear problems.⁵⁹ 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 Eq. (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 W_d(x), provided that it vanishes near all the boundaries. In particular, solutions which satisfy Eq. (6) can be computed using Eq. (11) combined with the relaxation process in which $ρ→1$ (or d → ∞). This idea is similar to the damping filter method.⁵⁸ 

Whether weighting is used or not, the Newton-Krylov solver generates the Floquet multipliers Λ_i and Floquet modes e_i of the computed solution, i.e., eigenvalues and eigenmodes of the full-space Jacobian $J=∂u(UTu)$⁠,

essentially “for free.” Indeed, the spectrum of the Krylov-subspace Jacobian H_k yields a good approximation to the leading eigenvalues and eigenmodes of A, while $J=g−1(A+1)$⁠. Of particular importance are the unstable modes (which correspond to |Λ_i| > 1) and the Goldstone modes (which correspond to |Λ_i| = 1) that characterize the symmetries of the system.

In this section, we present unstable spiral wave solutions of the modified Karma model computed using the weighted Newton-Krylov method described in Sec. III and investigate their properties. In particular, the spiral wave solution shown in Fig. 3(a) was found using continuation of a stable single-spiral solution with an approximately centered tip [x, y] = [74.43, 96.04], achieved by decreasing β (which corresponds to increasing the restitution parameter R). It has a period T = 50.8273 and wavelength λ = 74, or 1.94 cm in dimensional units. As Fig. 3(b) shows, the solution has one Goldstone mode and three unstable modes.

The Goldstone mode presented in Fig. 3(f) corresponds to the temporal derivative of the solution ∂_tu, as expected for a periodic orbit. The real unstable mode (Fig. 3(g)) corresponds to an almost rigid shift of the spiral wave in the y direction and can be identified with a frustrated translational Goldstone mode ∂_yu (its amplitude varies in space). The stable eigenvalue Λ = 0.6882 also corresponds to a frustrated Goldstone mode (Fig. 3(e)), which can be identified as linear combination of translations in the x and y directions, $n̂·∇u$⁠, where $n̂$ is the direction of translation. The two complex conjugate unstable modes (Figs. 3(c) and 3(d)) 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 R > 1. Since the Goldstone modes associated with spatial translations are frustrated for s = 32 (continuous translational symmetry is broken), no additional constraints beyond orthogonality with respect to ∂_tu are needed. The residual is minimized for h = 0, 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²³ for the modified version of the 3-variable Fenton-Karma (3V-FK) model⁶⁰ (they also replaced the Heaviside step functions with smoothed versions Θ_s(u)). The unstable spiral waves of 3V-FK (described by relative equilibria) were found to possess two near-Goldstone modes corresponding to spatial translations with Λ ≈ 1 and one Goldstone mode corresponding to temporal translation (or spatial rotation) with Λ = 1. 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 ∂_xu, ∂_yu, and ∂_tu, which become Goldstone modes in the presence of global continuous symmetries. For the Goldstone modes, we should have $(J−1)ei=0$ (up to the level of precision defined by ε_tol). We checked this by evolving the group tangents using the linearization of Eq. (1) about the periodic solution shown in Fig. 3(a). As Fig. 4 illustrates, the temporal tangent is a Goldstone mode: $|(J−1)∂tu|=O(10−11)$ while the spatial tangents are not. Both $(J−1)∂yu$ and $(J−1)∂xu$ achieve their O(10⁻²) maxima on the boundary of Ω (Figs. 4(d) and 4(e)), 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 ∂_yu and ∂_xu behave like Goldstone modes everywhere except near the core of the spiral wave (Figs. 4(g) and 4(h)), suggesting that translational symmetry is also broken locally. We investigate this local mechanism next.

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, Δx = Δy = 1. The spiral wave solution is characterized by several scales: the wavelength λ = 74, the length scale $ℓu=D11≈2$ of the fast (voltage) variable, which defines the width of the sharp leading front of the excitation wave, and the length scale l_v over which the dynamics of the slow variable v switches from excitable to refractory. The latter length scale is controlled by the term βΘ_s(u − 1) in Eq. (2) and, to leading order in ϵ, is given by l_v = 2ℓ_u/(sβ). Setting l_v = 1 gives s = 2ℓ_u/β ≈ 3. Hence, for l_v ≲ 1 (s ≳ 3) the solution is spatially ill-resolved, and the continuous translational symmetry is broken near the tip, pinning the spiral. On the other hand, for l_v ≳ 1 (s ≲ 3) 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 s = 1.257. The spiral wave solution at this low value of s is similar to that shown in Fig. 3(a), but corresponds to a generalized relative periodic orbit with period T = 54.7447 and wavelength λ = 78. The displacement of the wave over one period is small, but non-zero (⁠$|h|=O(10−9λ)$⁠), and its direction depends on both the position and the phase of the spiral wave. The corresponding spectrum (Fig. 5(a)) contains two complex conjugate unstable eigenvalues which correspond to the alternans instability (the corresponding modes are similar to those shown in Figs. 3(e) and 3(f)) and three eigenvalues with $|Λ−1|=O(10−6)$⁠. The corresponding Goldstone modes, predictably, coincide with the three group tangents ∂_tu (Fig. 5(d)), ∂_yu (Fig. 5(c)), and ∂_xu (Fig. 5(b)). 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-s regime where local Euclidean symmetry is preserved and the large-s regime where it is broken is continuous. The s-dependence of the leading eigenvalues is shown in Fig. 6(a). As the value of s is increased, two of the three unit eigenvalues split off around s = 3 and separate along the real axis. For this particular solution, Goldstone mode ∂_yu becomes unstable (and frustrated), while the Goldstone mode ∂_xu becomes stable (and frustrated) at large s. The exact symmetry of the problem with respect to rotations by $ϕ=π/2$ 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 x and y 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 1.2 ≲ s ≲ 6, 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 s.

The variation of h, which defines the net spatial displacement of the wave over one period of rotation, with s is shown in Fig. 6(b). Consistent with our dimensional analysis, we find that h vanishes for s ≳ 3, but does not vanish for s ≲ 3. Even though the spirals drift for small s, 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 s, as Fig. 6(c) illustrates. The deviation from the limiting value T₀, which corresponds to s → ∞, was found to be well-approximated by a power law, T − T₀ ∝ s^−4∕3. The overall variation, however, was not large, with the period being just 12% larger at s = 1 than at s = 32.

To further verify our dimensional analysis, we also computed the single-spiral solution on progressively finer grids with fixed stiffness parameter s = 32 and physical domain size L = 192 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 $|Λ−1|=O(10−4)$ when Δx = Δy = 0.16, whereas the dimensional analysis predicts a similar critical length scale l_v = 0.09 for this value of s.

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 ∂_tu = 0) has coordinates [x, y], then x and y 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 x) 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-s regime local continuous symmetry is broken, so we only find solutions shifted by an integer multiple of the grid spacing Δx. Regardless of the direction (towards the left or right wall), the spiral wave solution with y ≈ L/2 shown in Fig. 3(a) could only be continued for $|x−L/2|≲0.70λ$⁠, i.e., no closer than 0.60λ to either wall (cf. Fig. 7(c)). We will denote this solution branch u₊. When x was decreased below about 0.60λ, Newton-Krylov solver converged to a nearby spiral wave with y ≈ L/2 − Δy/2. This new spiral wave solution can be continued until x ≈ 0.33λ. We will denote this solution branch u₀. If continued in the opposite direction, u₀ can be extended symmetrically to x ≈ L − 0.33λ.

An example of a spiral wave corresponding to the branch u₀ with the tip near [x, y] = [96.5, 95.5] along with its spectrum is shown in Fig. 8. Its shape is essentially indistinguishable from that of u₊ (cf. Fig. 3(a)). Its spectrum is also similar to that of u₊ (cf. Fig. 3(b)), except for the eigenvalues with the positive real part: unlike u₊ which has a pair of eigenvalues, one unstable and one stable, on the real axis, u₀ has two stable complex conjugate eigenvalues Λ_± = 0.4662 ± 0.1735i. The real and imaginary parts of the corresponding modes are shown in Figs. 8(e) and 8(f) 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 parts of the complex conjugate pair of unstable modes are shown in Figs. 8(c) and 8(d) and correspond to the alternans instability.

Counting u₀, u₊, and $u−=Rπ/2u+$⁠, there are at least three distinct pinned spiral wave solutions in the large-s limit. Their tips are located, modulo Δx, approximately at [0.5, 0.5] for u₀, [0.5, 0] for u₊, and [0, 0.5] for u_−. 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 Δx, there are O(N²) “copies” of each of these three solutions.

Most of these “copies” are nearly identical. Consider, for instance, the solution u₀. Its leading eigenvalues (cf. Fig. 7(a)) are effectively independent of the position of the spiral wave and only begin to vary noticeably for x ≲ 0.4λ. Similarly, the period of this solution (cf. Fig. 7(b)) is essentially independent of the position of the tip over almost the entire range of x. The deviation of the period from the reference value T₀ at x ≈ L/2 decreases exponentially fast with x: $|T−T0|∝ exp(−x/ℓc)$⁠, where ℓ_c ≈ λ/24 is the characteristic length scale that describes the interaction of this spiral with the boundary.

In the small-s regime, the continuous translational symmetries are recovered, so there is a continuum of drifting spiral wave solutions parametrized by the coordinates x = [x, y] of the tip. Unlike the large-s limit, there is only one type of solution. Its leading eigenvalues Λ, the spatial shift h over one period, and the period T as a function of x (with y = L/2) are shown in Fig. 9. Just like in the large-s limit, we find all the basic properties of the spiral wave solution to be effectively independent of the position of the tip, provided x ≳ 0.5λ. The differences between distinct spiral wave solutions are exponentially small. For instance, Figs. 9(b) and 9(c) show that $|h|∝ exp(−x/ℓc)$ and $|T−T0|∝ exp(−x/ℓc)$⁠, where T₀ is the period of the centered solution and now ℓ_c ≈ λ/20. The largest differences (⁠$|h|≈0.05λ$ and $|T−T0|≈0.02T0$⁠) correspond to the distance of the closest approach x ≈ 0.14λ = O(ℓ_c).

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 s or continuously for small s) without changing any of its properties up to some level of resolution ε_res. 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 ε_tol. 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-s limit. As Fig. 9(a) illustrates, the eigenvalue associated with the x-translation deviates from unity for x ≲ 0.6λ. The vertical boundary does not break the y-translation symmetry for spiral waves (and hence does not cause a significant deviation of the corresponding eigenvalue from unity), unless they drift in the x direction. This symmetry is eventually broken for x ≲ 0.3λ when displacement h_x becomes significant. Comparison with Fig. 9(a) shows that for the corresponding eigenvalue $|Λ−1|=O(|hx|/λ$⁠).

It is worth noting that h_x > 0 for x ≲ 0.27λ, and h_x < 0 for 0.27λ ≲ x ≲ 0.40λ, 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^61,62 investigated the interaction with the boundaries for resonantly driven (and hence drifting) stable spiral waves in the Barkley model⁵ of excitable media. They also found that the interaction is repulsive at close range and showed that the interaction length scale ℓ_c is determined by the spatial extent of the response functions (adjoints of the Goldstone modes), which are localized to the core region.^50,51,63 We have confirmed by direct calculation that both the spatial and temporal response functions in the Karma model (not shown) decay as $e−ϱ/ℓc$⁠, where ϱ is the distance to the tip.

Finally, the period of rotation shows opposite trends for pinned (large s) and drifting (small s) spirals. As a pinned spiral approaches a boundary, it rotates more slowly (T increases), whereas a drifting spiral rotates more quickly (T 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.

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 s or discrete for large s) 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 Δx and varying physical size L. 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 u₁ and u₂ computed on domains of size L₁ = 448 and L₂ = 432, respectively. The difference δu = u₁ − u₂ inside Ω₂ is found to be concentrated in a narrow boundary layer of width O(ℓ_c) (cf. Fig. 10(b)), where ℓ_c ≈ λ/20 for s = 1.257, as we determined previously. Furthermore, Fig. 10(a) shows that the difference becomes exponentially small away from the boundaries, $|δu|∝ exp(−r/ℓc)$⁠, where r 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 s = 32, the spiral waves are pinned, and we find a minimal domain size, L₀ = 17 ≈ 0.23λ, approximately 0.45 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 L₀/2 = 0.115λ ≈ 2ℓ_c ≈ λ/12. Since ℓ_c determines both the radius of the spiral wave core and the width of the boundary layer, L = L₀ corresponds to the collision of the spiral core with the boundary layer. If this criterion applies more generally, we should expect to find spiral wave solutions only on domains of size L ≳ 4ℓ_c.

Pinned spiral waves remain unstable in the entire range of system sizes. The dominant eigenvalues of the solution u₀ are shown in Fig. 11(a). As L 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 L₀ ≤ L ≤ 0.36λ, there are between three and five unstable modes. As L → L₀, the dominant eigenvalues quickly grow in magnitude, approaching Λ_± = 2.0 ± 2.77i (these are outside of the range of Fig. 11(a)).

Interestingly, the character of the instability changes for L < L_b ≈ λ. Consider, for instance, the solution at L = 48 ≈ 0.65λ shown in Fig. 12(a). The magnitude of its unstable eigenvalues (cf. Fig. 12(b)) is comparable to that in much larger systems, but the corresponding eigenmodes (Figs. 12(c) and 12(d)) correspond to the meandering instability, rather than alternans. Initial conditions close to this unstable solution produce spiral waves which persist for up to 10³ rotations without breaking up, which is consistent with numerical results of Karma⁶ for domains smaller than 2.1 cm (which corresponds to L_b = 80 = 1.08λ 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 L_b that supports spiral wave breakup via alternans defines another important dynamical length scale.

The period of the solution approaches the asymptotic value T₀ exponentially fast as L increases (cf. Fig. 11(b)). For L ≲ 3λ, we find $|T−T0|∝ exp(−L/2ℓc)$⁠, which is consistent with the scaling we found previously. Indeed, for a domain-centered spiral wave x = L/2, or L = 2x, which means that scaling with x and L follows from the same general scaling law. For L ≳ λ, the period can be considered essentially independent of L.

Most of the results for small s 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 L, $|h|∝ exp(−L/2ℓc)$ (cf. Fig. 14(b)). As the domain is truncated, the eigenvalues associated with Goldstone modes begin to deviate from unity (cf. Fig. 14(a)). Eigenvalues corresponding to translational modes deviate first, at L ≈ 1.3λ. This agrees well with the distance x ≈ 0.6λ to the boundaries at which translational symmetry is broken, as we found in Sec. IV B. Somewhat unexpectedly, for even smaller domains (L ≲ 0.6λ), the eigenvalue associated with temporal translation also deviates from unity. The critical domain size matches the distance of the tip to the boundary (x ≈ 0.3λ) beyond which the spiral solution loses translational symmetry, as noted in Sec. IV B. 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 L ≲ 0.6λ, 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 u(0) and u(T) become noticeably different.

The decrease in the temporal period at small separations between the tip of the spiral and the boundaries observed in Figs. 9(c), 11(b), and 14(c) is due to the negative curvature of the wave front near the boundary due to the no-flux boundary condition. It is well-known⁶⁴ that the speed of propagation of an excitation wave is related to the curvature κ of its front by a linear relationship, with convex (κ > 0) wave fronts traveling more slowly than flat (κ = 0) ones, and flat wave fronts traveling more slowly than concave (κ < 0) ones.^65,66 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 ℓ_c.

The exponential dependence of the shift magnitude $|h|$ and the period T of the drifting spiral waves found for the modified Karma model are in contradiction with the analytical results obtained by Aranson et al.⁶⁷ for a similar model in the ν → 0 limit, which predict a super-exponential dependence of the drift speed $|h|/T$ and rotational frequency ω = 2π/T of spiral waves on the domain size L. Davydov and Zykov⁶⁸ predict the frequency of spirals in a generic reaction-diffusion model with ν = 0 to vary as the inverse of the domain size (for circular domains of radius comparable to λ). Hartmann et al.⁶⁹ 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 are in excellent agreement with the exponential dependence.

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 O(ℓ_c), 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 O(ℓ_c).

Global Euclidean symmetry can, of course, be gradually restored by increasing the size L 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 L increases. The difference between solutions with different L is again found to be confined to the boundary layer of width O(ℓ_c). Similarly, the properties of solutions approach the asymptotic values exponentially fast as L increases, with the same length scale ℓ_c. This length scale (ℓ_c) was found to control the effect of the boundaries rather generally. Effectively, ℓ_c 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 $exp(−r/ℓc)$⁠, with the distance r 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, L ≥ L₀, with the smallest domain size corresponding to the strongly interacting regime L₀ ≈ 4ℓ_c 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 s considered here, the minimal spacing between spiral cores in the state of fibrillation.⁴⁸ As L varies between L₀ and λ/2, 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.⁷⁰ 

We have established that, for L < L_b meander is the leading instability mechanism, so spirals of size L₀ < L < L_b 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 L > L_b 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 L₀). 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 L₀ to L_b.

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.