Unstable nonchaotic solutions embedded in the chaotic attractor can provide significant new insight into chaotic dynamics of both low- and high-dimensional systems. In particular, in turbulent fluid flows, such unstable solutions are referred to as exact coherent structures (ECS) and play an important role in both initiating and sustaining turbulence. The nature of ECS and their role in organizing spatiotemporally chaotic dynamics, however, is reasonably well understood only for systems on relatively small spatial domains lacking continuous Euclidean symmetries. Construction of ECS on large domains and in the presence of continuous translational and/or rotational symmetries remains a challenge. This is especially true for models of excitable media which display spiral turbulence and for which the standard approach to computing ECS completely breaks down. This paper uses the Karma model of cardiac tissue to illustrate a potential approach that could allow computing a new class of ECS on large domains of arbitrary shape by decomposing them into a patchwork of solutions on smaller domains, or tiles, which retain Euclidean symmetries locally.

Spiral turbulence in excitable systems, which features multiple interacting spiral waves that repeatedly break up and merge, is of substantial practical interest because of its relation to life-threatening cardiac arrhythmias such as ventricular fibrillation. While mathematical models of cardiac tissue vary in complexity, even the simplest, such as the Karma system, are too high-dimensional to yield much insight, limiting our understanding of the dynamical mechanisms that sustain spiral turbulence. A growing body of recent work in fluid turbulence supports the conjecture that high-dimensional, spatiotemporal chaos can be understood as a series of transitions between unstable recurrent patterns, also known as exact coherent structures, embedded in the chaotic set on which these dynamics take place. If successful, a similar approach could allow construction of a dynamical description of spiral turbulence in cardiac tissue in terms of models with dramatically reduced dimensionality. Such a reduced-order description may significantly advance our understanding of spiral turbulence and pave the way towards more effective defibrillation techniques. However, the local Euclidean symmetries, inherited by spiral waves from the respective global symmetries of the underlying evolution equations, significantly complicate the job of finding any exact coherent structures. The objective of this paper is to develop a formalism for local symmetry reduction that would enable efficient computation of these structures.

## I. INTRODUCTION

Cardiac arrhythmias, such as atrial and ventricular fibrillation, are characterized by spatially complex, high-dimensional dynamics that are generated as multiple excitation waves propagate through cardiac tissue, merging and breaking up. Despite substantial advances in computing power and the development of detailed ionic models of cardiac cells, quantitatively accurate direct numerical simulation (DNS) of cardiac tissue remains computationally expensive and provides limited dynamical insight into mechanisms that initiate and sustain the spatiotemporal chaos that underpins these arrhythmias.

Substantial progress in understanding some types of spatiotemporal chaos has been made over the past two decades using an idea that is now over a century old. In developing celestial mechanics, Poincaré^{1} realized that unstable equilibria and periodic orbits provide a skeletal structure which organizes chaotic dynamics. His idea was later developed in the context of quantum chaos by Gutzwiller^{2} and subsequently applied to high-dimensional chaos generated by nonlinear partial differential equations (PDEs) such as the Kuramoto-Sivashinski equation^{3,4} and Ginzburg-Landau equation.^{5}

Although in one spatial dimension unstable periodic solutions could be computed using brute-force Newton iteration, this numerical approach becomes intractable for two- and three-dimensional PDEs whose discretizations routinely involve millions of degrees of freedom. In this case both nonchaotic solutions and their spectra can be computed efficiently^{6} using a combination of Newton descent, Krylov subspace/GMRES solution of the Newton equations, and “trust-region” heuristic for the magnitude of the Newton steps. Newton-Krylov methods facilitated recent studies of various fluid flows at intermediate Reynolds numbers,^{7–10} which produced valuable new insight into the mechanisms that generate and sustain fluid turbulence—arguably the most challenging unsolved problems of classical physics. Although periodic orbit theory has never been used to analyze spatiotemporally chaotic dynamics in excitable systems, its success in uncovering the mysteries of fluid turbulence gives us hope that it can also generate new insights into the problem of fibrillation and thereby help develop new and improved methods of defibrillation.^{11–13}

Despite the progress that has already been made in using unstable nonchaotic solutions to understand chaotic dynamics, many open problems remain. In particular, it is not always clear what types of unstable nonchaotic solutions play a dominant role in spatiotemporal chaos. For spatially extended systems, nonchaotic solutions are characterized not only by their temporal properties (e.g., equilibria, periodic orbits), but also by their spatial structure. In fact, the spatial structure received far more attention in the studies of fluid flows, which provided the motivation for the recent development of periodic orbit theory. As a result, nonchaotic solutions of nonlinear PDEs embedded in the chaotic set on which the dynamics take place have become known as exact coherent structures (ECS), reflecting their connection with the empirically observed coherent structures in fluid flows.

Spatially extended systems often respect continuous Euclidean symmetries, which complicate the dynamical description based on periodic orbit theory. In particular, continuous symmetries give rise to several other dynamically relevant classes of nonchaotic unstable solutions, such as relative equilibria and relative periodic orbits, which reduce to equilibria and time-periodic orbits in a co-moving reference frame. Notably, the numerical methods for computing such solutions in co-moving frames were developed in the context of excitable/oscillatory systems such as the Barkley model,^{14,15} FitzHugh-Nagumo, $\lambda \u2212\omega $, and complex Ginzburg-Landau model.^{16} In two (or three) dimensions, however, rotational symmetry requires that the computation be performed on a circular (or cylindrical) domain, using a polar grid, which severely limits the usefulness of this approach. We recently showed^{17} that, for systems with local Euclidean symmetries, relative equilibria and relative periodic orbits could be computed on domains of arbitrary shape, using arbitrary grids, with the help of a weighted Newton-Krylov method.

In this work, we use a prototypical model of cardiac tissue to determine which types of nonchaotic solutions form a skeleton for the chaotic set on which spatiotemporally chaotic dynamics that underlie fibrillation take place. Although numerous nonchaotic solutions involving few unstable spirals with well-separated cores were computed using the weighted Newton-Krylov method, these solutions were found to be located away from the chaotic set. Attempts to find unstable solutions embedded in the chaotic set, which feature multiple spirals with closely-spaced cores, were unsuccessful. The analysis of converged few-spiral solutions suggests that this failure is due to the local Euclidean symmetries which cause slow relative drift and/or rotation of the spiral cores. If this is indeed the case more generally, it may be possible to decompose multi-spiral solutions into tiles,^{18,19} each of which contains a single spiral described by either a periodic or a relative periodic solution. The ability to reduce local symmetries via such a decomposition will provide an important first step towards constructing a reduced-order dynamical description of spiral turbulence.

This paper is organized as follows. In Sec. II, we describe the modified Karma model of cardiac tissue used in this study. In Sec. III, we discuss the relationship between Euclidean symmetries and different classes of nonchaotic solutions. Section IV describes our failed attempts to find dynamically relevant nonchaotic solutions embedded in the chaotic set. In Sec. V, we explore the effects of local symmetries using artificially constructed unstable nonchaotic solutions with multiple spirals. The construction of tiles and their use for symmetry reduction are discussed in Sec. VI. Our conclusions and directions of further research are presented in Sec. VII.

## II. THE KARMA MODEL

Mathematical models of cardiac tissue vary greatly in complexity. The most detailed models of cellular kinetics include tens of ordinary differential equations describing the evolution of intracellular concentrations of various ions and the state of ionic channels as well as the transmembrane potential (i.e., the difference between the intra- and the extracellular electric potential). Monodomain tissue models also include the electrical coupling between the cells. Bidomain models go one step further, accounting for the spatial and temporal variation of both the intracellular and the extracellular potential.

The simplest monodomain models can have as few as two local variables that describe the transmembrane potential and the internal state of a cardiac cell (cardiomyocyte). Despite this simplification, they can produce dynamics that are as rich and complicated as those of the bidomain models including a detailed description of the cellular kinetics. Most monodomain models are in the class of reaction-diffusion systems described by coupled nonlinear PDEs of the form

where the vector field $w=(u,v)$ incorporates the trans-membrane voltage *u* and gating variable(s) v, *D* is a diagonal matrix of diffusion coefficients, *τ _{u}* is the characteristic time scale, and the nonlinear function $f(w)$ describes the cellular kinetics.

Although detailed multi-variable models may provide an accurate description of single cell dynamics, they often fail to accurately describe the dynamics of cardiac tissue. In addition, they require substantial computational resources to obtain numerical solutions with appropriate spatial and temporal resolution. Furthermore, when compared to some algebraically simpler, lower-dimensional models, the added complexity provides little additional insight into the *dynamical* mechanisms responsible for generating and maintaining complex arrhythmic behaviors such as fibrillation.

In contrast, the two-variable Karma model^{20} provides a greatly simplified description of cardiomyocyte excitation dynamics while still reproducing an essential alternans instability. This instability is believed to be responsible for the breakup of a single-spiral solution associated with the arrhythmic state of tachycardia and, ultimately, the transition to fibrillation in two or three dimensions, as illustrated in Fig. 1.

We use a slightly modified form of the Karma model for the cellular kinetics

where $\Theta s(u)=[1+tanh(su)]/2$. The parameter *ϵ* describes the ratio of excitation and relaxation time scales, $u*$ is a phenomenologically chosen voltage scale, and $\beta =1/(1\u2212exp(\u2212R))$ controls the restitution and susceptibility of traveling excitation waves to alternans. Following the original study,^{20} we choose the parameters $\tau u=2.5$ ms, $u*=1.5415,\u2009\u03f5=0.01,\u2009Du=0.11$ mm^{2}/ms, and *R* = 1.273 such that spiral wave solutions break up as a result of the alternans instability.

The modifications (the addition of weak diffusion in the gating variable, replacing the Heaviside step function $\Theta (u)$ with its smoothed version $\Theta s(u)$, and the additional term $\Theta s(v\u22121)(v\u22121)$) serve to smooth the unphysical discontinuities and singularities of the original model without noticeably changing its dynamics. The stiffness parameter *s* (in this study we set *s* = 32) controls the switching of the gating variable as well as the shortest length and time scales of the spiral wave solutions.^{17} Existing models of cardiac tissue, including the original Karma model, ignore diffusion in describing the dynamics of the gating variable(s). However, it is well-established, although perhaps not well-known, that all relevant ions and small secondary messenger molecules pass through the gap junctions between cardiac cells.^{21,22} Hence, weak diffusion (we set $Dv=5.5\xd710\u22123$ mm^{2}/ms) of the gating variable is physiologically justified.

Our focus here is on spiral turbulence in two dimensions, which can be thought of as a model of atrial fibrillation. Without loss of generality, we can assume that the computational domain Ω is rectangular. The PDEs (1) are therefore discretized on a rectangular uniform two-dimensional grid, with *N* grid points in the *x*-direction and *M* grid points in the *y*-direction. The use of a rectangular grid is consistent with the structure of cardiac tissue, whose properties are generally anisotropic (i.e., different in the direction along, and transverse to, the muscle fibers). We use a grid spacing $\Delta x=\Delta y=262$ *μ*m which corresponds to the typical size of cardiomyocytes and use the physiologically relevant no-flux boundary conditions $n\xb7\u2207w=0$, unless otherwise specified.

The equations are solved using a fourth-order Runge-Kutta time integrator and a nine-point finite difference approximation for the Laplacian that minimizes the effect of discretization on the rotational symmetry of the evolution equations.^{17} The time step used in the simulations was $\Delta t=0.01$ ms, except for the states shown in Fig. 4 where $\Delta t=0.1$ ms.

## III. UNSTABLE NONCHAOTIC SOLUTIONS AND EUCLIDEAN SYMMETRIES

On an unbounded two-dimensional domain $\Omega =\mathbb{R}2$, reaction-diffusion equations are equivariant under the set of transformations, or actions, *g* that form a group $G=G\u2297E+(1)$. In particular, *spatial* transformations from the Euclidean subgroup $G=E(2)$ include continuous translations, continuous rotations and discrete reflections. The subgroup $E+(1)$ reflects the fact that the dynamics are invariant under *temporal* translations, but are not time-reversible.

Spatially uniform solutions **w** such that $f(w)=0$ (e.g., the stable resting state $w\u2248(0,0)$), invariant with respect to the entire group *G*, are examples of the simplest type of nonchaotic solutions–equilibria. However, typical nonchaotic solutions respect only some of the transformations in *G*, lowering effective symmetry to a subgroup of *G*. In particular, a plane wave solution traveling with velocity **c** is a relative equilibrium generated by a continuous translational symmetry which satisfies $\u2202tw\u2212c\xb7\u2207w=0$. In this case $\u2202t$ is the generator of translations in time and $\u2207$ is the generator of translations in space. Another relative equilibrium generated by continuous rotational symmetry, which satisfies $\u2202tw\u2212\omega \u2202\theta w=0$, corresponds to a spiral wave. Here *θ* is the polar angle in the (*x*, *y*) plane, *ω* is the rotational frequency, and $\u2202\theta $ is the generator of rotations. These relative equilibria reduce to (nonuniform) equilibria in a reference frame translating with velocity **c** or rotating with angular velocity *ω*, respectively. For these special solutions, evolution (i.e., time translation) is equivalent to a spatial translation or rotation.

The introduction of physical boundaries has a markedly different effect on these solutions. For no-flux boundary conditions, the effect of the boundaries is determined primarily by the shape (and size) of the domain Ω on which the evolution equations are defined. Some nonchaotic solutions, such as the uniform equilibria or spiral waves, survive unchanged (or almost unchanged) when computed on domains of very different shapes and sizes, while others, such as plane waves, cease to exist except as a transient. The effect of the boundaries on spiral waves is of the most interest to us, since spiral waves play an important role in sustaining spatiotemporally chaotic dynamics underlying complex cardiac arrhythmias.

Although it is easier to compute single-spiral waves described by relative equilibria on circular domains, resulting solutions are not representative of any cardiac rhythms. For instance, a spiral wave describing tachycardia would not correspond to a relative equilibrium because atria are not circular. Similarly, the atria cannot be tiled with circles, so multi-spiral patterns arising during atrial fibrillation cannot be decomposed into relative equilibria, even on arbitrarily short time scales. Hence, it is crucial to understand the properties of spiral wave solutions on domains of arbitrary shape. Rectangular domains are more relevant, since they break all global Euclidean symmetries and, at the same time, can be used to tile a surface. Fig. 1(a) shows a sample unstable periodic solution computed on a square domain of side *L* = 50.3 mm with no-flux boundary conditions using the weighted Newton-Krylov solver.^{17} It describes a pinned spiral wave with a period of *T* = 125.98 ms and wavelength $\lambda =19.4$ mm (for reference, the left human atrium has a typical diameter of 30–40 mm).

Since this is not a relative equilibrium, the actions of time translation and rotation become distinct, $\u2202tw\u2260\omega \u2202\theta w$, where $\omega =2\pi /T$. However, the effect of the boundaries on the spatial structure of the spiral wave as well as its temporal properties (e.g., its period *T*) becomes exponentially weak for sufficiently large *L*. The solutions computed on domains of different size become indistinguishable to numerical precision in the overlap region, with the difference becoming significant only near the boundary.^{17} In particular, the periodic solution shown in Fig. 1(a) is almost identical to a relative equilibrium in the interior of the domain. Hence, although the boundary conditions break the global rotational symmetry of the spiral wave, local rotational symmetry is preserved away from the boundaries, where $\u2202tw\u2248\omega \u2202\theta w$ to numerical accuracy.

For the parameters chosen in this study, single-spiral solutions are unstable with respect to the alternans instability. If allowed to evolve for a sufficiently long time, an initial condition that is arbitrarily close to the spiral wave shown in Fig. 1(a) will deviate from that unstable solution exponentially fast, developing modulation in the width of the excitation wave (alternation of the action potential duration), as illustrated in Fig. 1(b). When the amplitude of the modulation becomes sufficiently large, the tissue fails to recover between two successive wave fronts, leading to conduction block and a breakup of the spiral wave (cf. Fig. 1(c)). After a sequence of breakups, a state of spiral turbulence featuring multiple interacting spiral waves is established (cf. Fig. 1(d)).

In the process of computing a solution, the Newton-Krylov solver also generates its stability spectrum. The leading Floquet modes (or simply eigenmodes) of the single-spiral solution depicted in Fig. 1(a) are shown in Fig. 2. The pair of unstable eigenmodes on the left, corresponding to a complex conjugate pair of Floquet multipliers (or eigenvalues) outside the unit circle, are associated with the variation in the width of the excitation wave (i.e., alternans instability). A single unstable eigenmode on the right looks like a frustrated Goldstone mode $\u2202yw$; it is associated with a weakly broken translational symmetry and corresponds to an almost rigid vertical shift of the entire spiral wave. The Goldstone mode $\u2202tw$ has a unit eigenvalue Λ = 1 and is associated with a continuous symmetry (time-translation).

For a spiral wave solution of the PDEs (1) on an unbounded domain, we should expect to find three Goldstone modes corresponding to each of the three continuous symmetries of *E*(2): rotation and translations in the *x* and *y* directions. The continuous time-translation symmetry does not produce an additional Goldstone mode because for spiral waves described by relative equilibria $\u2202tw=\omega \u2202\theta w$. On finite, but sufficiently large, domains the three continuous spatial symmetries become local, rather than global, but should still generate three unit eigenvalues, even after (1) is discretized, provided the grid is sufficiently fine. Indeed, this was found to be the case^{17} for small values of the stiffness parameter ($s\u2009\u2272\u20093$). As is typical of solutions in the presence of continuous translational symmetries, in this limit the spiral waves drift, and hence are described by relative periodic solutions (they become periodic in a reference frame moving with the drift velocity).

However, for the values of parameters considered here, there is also a local source of continuous symmetry breaking – the discreteness associated with the computational grid (or the cellular structure of the tissue). The modified Karma model possesses an internal length scale $\u2113v=2Du\tau u/(\beta s)$. For *s* = 32, this gives $\u2113v\u223c24$ *μ*m, i.e., much smaller than the grid spacing $\Delta x=262$ *μ*m. Hence the computational grid is too coarse to completely resolve the spatial structure of the spiral wave, effectively reducing the continuous translational symmetries to discrete translations, $x\u2192x+\Delta x$ and $y\u2192y+\Delta y$, such that $\u2202xw$ and $\u2202yw$ no longer correspond to eigenmodes.^{17}

Our numerical results are consistent with this conclusion. On square domains of sufficiently large size, for $s\u2009\u2273\u20093$, there is a large, but finite, number of unstable spiral waves described by periodic solutions with essentially identical periods and stability spectra and differing by a discrete shift in the *x* and/or *y* direction. Fig. 1(a) shows just one of those solutions. Each one of these spiral waves possesses a Goldstone mode $\u2202tw$ which coincides with $\u2202\theta w$ (to numerical precision) away from the boundaries, i.e., continuous rotational symmetry is preserved locally.

## IV. EXACT COHERENT STRUCTURES

Stable and unstable uniform states, plane and spiral excitation waves discussed previously illustrate the types of solutions (equilibria, periodic orbits, relative equilibria, and relative periodic orbits) relevant to cardiac tissue dynamics, including fibrillation. However, none of these solutions are embedded in the chaotic set on which arrhythmic dynamics underlying fibrillation take place, i.e., they do not correspond to ECS. For instance, while locally the excitation waves almost always take the shape of small spirals during fibrillation, they never organize globally into a single large spiral wave.

In order to find ECS, i.e., global unstable nonchaotic solutions embedded in the chaotic set, we used the method of close returns^{23} which has been used successfully in the context of fluid turbulence.^{24} The procedure involves finding near-recurrences in the *chaotic* solutions, which become initial conditions that are subsequently refined into exact *nonchaotic* solutions using the Newton-Krylov solver.^{17} In the presence of global symmetry, initial guesses for relative equilibria, periodic orbits, as well as relative periodic orbits can be found as the minima of the recurrence function

where $||\xb7||2$ denotes the 2-norm in $\mathbb{R}2NM$ and the group $G$ describes the symmetries in the presence of boundaries. However, we found that the minima of (3) are always achieved for $g\u22481$, even for periodic boundary conditions, and in practice we can set $G={1}$.

It should be noted that setting $G={1}$ in (3) does not constrain the exact solutions to absolute equilibria and absolute periodic orbits. For instance, a slowly drifting or rotating solution described by a relative periodic orbit or relative equilibrium, such that $gw(t)=w(t\u2212T)$ with $g\u22481$, will generate a minimum of (3) with $\tau \u2248T$ when $G={1}$. Solutions that are characterized by fast global rotation or translation (e.g., spiral or plane waves), and for which *g* is significantly different from $1$, are characterized by a high degree of spatial coherence. The lack of global spatial coherence is a distinguishing feature of fibrillation, and therefore we should not expect to find any solutions that exhibit fast global rotation or drift.

A fragment of the recurrence plot for a chaotic solution computed on a square domain Ω with side *L* = 50.3 mm and periodic boundary conditions is shown in Fig. 3. The periodic boundary conditions lead to chaotic dynamics that are qualitatively identical to those in the presence of no-flux boundary conditions even on relatively long time scales. However, for periodic boundary conditions fibrillation persists indefinitely (at least, it does not disappear after 4.5 min $\u22482000T$), while for no-flux boundary conditions it can terminate spontaneously due to the spiral cores colliding with the boundaries and disappearing. Once a sufficiently low minimum (circled) of the recurrence function $E(t,\tau )$ is identified, the corresponding state $w(t\u2212\tau )$ is used as the initial guess for a solution with period close to *τ*.

Rather surprisingly, none of the initial guesses we tried converged to either relative equilibria, periodic orbits, or to relative periodic orbits. (We did not search for equilibria, since $||\u2202tw||$ never becomes small for chaotic solutions.) Figs. 4(a) and 4(b) show the voltage component *u* for two typical examples of the numerous multi-spiral states identified using the recurrence analysis. They nearly recur after $\tau =251.94$ ms and $\tau =251.83$ ms, respectively. These values of *τ* correspond to approximately double the temporal period of a single spiral shown in Fig. 1(a), $2T=251.96$ ms. Newton iterations stagnate at the values of relative residual $||w(\tau )\u2212w(0)||2/||w(0)||2$ equal to $9\xd710\u22123$ for the state shown in Fig. 4(a) and $10\u22122$ for the state shown in Fig. 4(b), compared to $O(10\u221213)$ for the converged single-spiral state shown in Fig. 1(a). The voltage components $u(t)\u2212u(t\u2212\tau )$ of the corresponding residual are shown in Figs. 4(c) and 4(d). They are spatially localized in the regions where $|\u2207u(t)|$ has the largest magnitude (cf. Figs. 4(e)–4(f)), or near the front and back of the excitation wave. In order to interpret these findings, next we consider multi-spiral solutions that are constructed artificially and do not lie on the chaotic set on which fibrillation takes place.

## V. MULTI-SPIRAL SOLUTIONS WITH LOCAL SYMMETRY

The weak effect of the boundaries on the essential dynamics of the spiral wave suggests that the local translational and rotational symmetries of single-spiral solutions should be inherited by multi-spiral solutions, provided the spacing between the spiral cores is sufficiently large. To verify this and to gain insight into the properties of multi-spiral states on finite domains, we explore the simplest artificially generated solutions, which are constructed by assembling pairs of converged single-spiral solutions.

### A. Local rotational symmetries

We first explore local rotational symmetry by constructing a set of initial states containing two spirals that are phase shifted with respect to each other. These initial states are prepared by placing a copy of a single, isolated spiral solution such as the one shown in Fig. 1(a) next to a second identical copy that has been integrated forward in time by a fraction of its temporal period. Both co-rotating and counter-rotating spiral states are prepared to account for possible effects of chirality. These initial states are then refined into exact solutions using the Newton-Krylov solver. Converged solutions for four phase shifts $\delta \theta ={0,\pi /2,\pi ,3\pi /2}$ are shown in Fig. 5 (counter-rotating) and Fig. 6 (co-rotating).

The leading eigenvalues for the co-rotating spirals with $(\delta \theta =\pi )$ shown in Fig. 6(c) are plotted in Fig. 7. This spectrum is representative of the spectra found for the other phase-shifted two-spiral solutions (both co- and counter-rotating) in the sense that all leading eigenmodes could be spatially localized to the region occupied by one or the other spiral. Each pair of localized eigenmodes corresponds to one of the eigenmodes shown in Fig. 2 for the single-spiral solution. The corresponding eigenvalues are degenerate (to numerical precision $O(10\u22126)$), indicating that the dynamics of the two spirals are independent. The pair of Goldstone modes associated with the two unit eigenvalues (Λ = 1) reflects the freedom of either spiral to undergo small phase shifts relative to the other without destroying the time-periodic nature of the solution globally. Just like for single spirals, in the interior of the domain the phase shifts are indistinguishable from rotation of one spiral relative to the other: rotation by angle $\delta \theta $ corresponds to a temporal shift $\delta t=(\delta \theta /2\pi )T$.

### B. Local translational symmetries

We explore local translational symmetries, which become discrete on our computational grid, by using two distinct sets of solutions. The first set is constructed by taking two single, isolated spiral solutions (in phase) and applying discrete shifts in the horizontal or vertical direction to one of the spirals by an integer multiple of the grid spacing before recombining them. Figure 8 shows a set of converged two-spiral solutions in which the core of the spiral on the right side of the domain has been shifted with respect to the spiral core on the left side of the domain. The shift for the spiral in Fig. 8(a) is by 30 grid points along the negative *y* direction (downward), while the shift for the spiral in Fig. 8(b) is 30 grid points in the negative *x* direction (leftward). While the convergence of these two solutions confirms that multi-spiral solutions respect locally discrete translational symmetries, the Goldstone modes associated with continuous translational symmetry remain absent in the spectrum for both solutions. Solutions with smaller shifts (for example, 5 grid points instead of 30) also converge and have similar spectral properties.

The second type of solutions used to explore local translational symmetry is shown in Fig. 9(a). It involves a three-spiral configuration in which a small spiral wave is formed between two larger spiral waves. Under time evolution, the small spiral undergoes a drift in the mostly vertical direction. This drift is caused by interaction with the neighboring spirals, as we explain below. Drifting spirals in bounded geometries are associated with relative periodic solutions which respect (local) continuous translational symmetries.^{17} The two big spirals do not drift, so the Newton-Krylov solver fails to converge onto a time-periodic (or relative periodic) solution, with the smallest relative residual magnitude equal to $4\xd710\u22123$. Figure 9(d) shows that the voltage component of the residual (with $\tau =125.98$ ms) is concentrated in the vicinity of the drifting spiral and (just like for the multi-spiral states shown in Fig. 4) is localized near the front and back of the excitation wave (cf. Fig. 9(b)), which is consistent both with the slow drift and the slow rotation of the small spiral.

This simple three-spiral solution is important since it is qualitatively representative of spiral interactions that take place during more complex chaotic behavior. It provides valuable insight into the dynamics and local symmetries of multi-spiral solutions such as those shown in Figs. 4(a) and 4(b). While each of the three individual spirals can be represented by either periodic or relative periodic solutions in the vicinity of its core, the solution on the entire domain is neither a periodic nor a relative periodic orbit. Similarly, one may expect multi-spiral ECS to be described by periodic and/or relative periodic solutions locally, but not globally.

## VI. STRATEGY FOR LOCAL SYMMETRY REDUCTION

A new class (or classes) of nonchaotic solutions are required to define a skeleton for the chaotic set associated with spiral turbulence. These nonchaotic solutions are associated with near-recurrences such as those shown in Figs. 4(a) and 4(b) that are found frequently, and hence play an important dynamical role, in spiral turbulence. In this section, we discuss how such nonchaotic solutions can be constructed starting with nearly-recurrent solutions that can be identified directly by exploring the natural measure on the chaotic set. As we have shown previously, spatial coherence in spiral turbulence is limited to areas (or tiles) associated with one spiral. Hence, we start by describing the process of decomposing the computational domain into individual tiles, each of which contains a single-spiral solution. We then discuss which boundary conditions these local single-spiral solutions should satisfy and how they can be assembled into a global nonchaotic solution.

### A. Tiling the domain

The concept of domain tiling was introduced by Bohr *et al.*^{18,19} to describe frozen spiral waves, sometimes referred to as vortex glasses, in the complex Ginzburg-Landau equation (CGLE)

where *A* is a complex field and *α* and *β* are control parameters. Each tile contains exactly one spiral core and the dynamics on each tile is controlled almost entirely by that core. The boundaries of individual tiles were identified with the ridges (or shocks) of the field $|A|$ which describes the amplitude of local oscillation. Fig. 10(a) shows the real part of a representative solution $A=\rho ei\varphi $ which can be used to identify individual spirals. In most of the domain, the phase $\varphi $ of the oscillation varies slowly in space, so according to the amplitude equation^{25}

the amplitude of oscillation is essentially constant, $\rho \u22481$, which corresponds to the middle cycle in the complex-*A* plane shown in Fig. 10(b). The spiral cores are associated with phase singularities and are characterized by small values of the amplitude ($\rho \u226a1$, small cycle). Similarly, the phase changes quickly at the boundaries of the tiles, where the amplitude increases ($\rho \u2009\u2273\u20091$, large cycle). Hence, the maxima and minima of $|A|$ shown in Fig. 10(c) can be used to identify, respectively, the spatial locations of the tile boundaries and spiral cores.

For excitable media characterized by strongly nonlinear oscillations, a different representation has to be used, since the phase and amplitude of the oscillations can be difficult to define, let alone compute. In this case the local amplitude of oscillation can be characterized instead by the area *I*(*x*, *y*) of the cycle *C* that is traced out by the solution in an appropriate state space. For CGLE the area in the complex-*A* plane is given by

with the result shown in Fig. 10(d). For frozen spirals $\varphi \u0307=2\pi /T$, so that $I(x,y)=\pi |A(x,y,t)|2$ and the cycle area representation is equivalent to the amplitude representation.

A similar approach can be used to identify the tile boundaries for the modified Karma model. The easiest way to characterize the amplitude of the strongly nonlinear oscillations is by computing the area *I*_{1} of the cycle *C* in the (*u*, *v*) plane,

Several representative cycles *C* computed for the co-rotating two-spiral solution depicted in Fig. 6(c) are shown in Figure 11(a) and the corresponding spatial distribution $I1(x,y)$ is shown in Figure 11(b). The two spirals are separated by a shock which corresponds to the local maximum of $I1(x,y)$. In addition, we also find shocks that form along the outer boundary, where no-flux boundary condition is imposed. Taken together, the shocks form a closed boundary for each of the spiral domains, defining the tiles on which the dynamics is controlled by one or the other core.

While computing cycle areas in the (*u*, *v*) plane is convenient in the models where all variables are accessible, in experiment this is rarely the case. Most typically only one variable is easily accessible experimentally (e.g., voltage or calcium, if voltage- or calcium-sensitive dye is used). In this case, cycle areas can also be computed using alternative planar representations of the cellular dynamics based solely on one variable, e.g., the voltage *u*. One possibility is to use a $(u,u\u0307)$ plane, with the cycle area defined as

Some representative cycles and the cycle area distribution are shown in Figs. 11(c) and 11(d), respectively.

The expression for *I*_{2}, however, involves a derivative of *u*. Since experimental measurements are typically noisy, derivatives obtained using finite-differencing of a time series can be very inaccurate. To reduce the influence of noise, a time-delay embedding can be used instead, with the cycles defined in the plane spanned by *u* evaluated at times *t* and $t\u2212\tau $ (the latter will be denoted $u\tau $ for brevity). The corresponding cycle area

can be computed without evaluating derivatives of any field. The choice of the time delay *τ* is not unique. We chose the value $\tau =13.5$ ms which corresponds to the first minimum of the mutual information function.^{26} The corresponding representative cycles in the $(u,u\tau )$ plane and the cycle area distribution are shown in Figs. 11(e) and 11(f), respectively. Comparison of Figs. 11(b), 11(d), and 11(f) shows that all three cycle area representations are consistent and accurately capture the shock line separating the two spirals, as well as the shock lines along the domain boundaries.

In the rest of this study, we use the cycle areas *I*_{1} computed in the (*u*, *v*) plane to identify tiles in the computational domain. In particular, Figs. 4(g) and 4(h) show the shocks that separate the tiles that form for the nearly-recurrent multi-spiral solutions, while Fig. 9(b) shows the shocks that separate the tiles in the three-spiral solution. For time-periodic solutions (e.g., single- or two-spiral solutions of the Karma model) the period *T* is well-defined and the cycles close perfectly. For non-periodic solutions such as those shown in Figs. 4 and 9, the integration is instead performed between crossings of a convenient Poincaré section (we used *u* = 1.5).

If the phase of the spiral solution is well-described by the Archimedian approximation, the tile boundaries can also be constructed analytically. Generally, the Archimedian approximation, and hence the analytic solutions for the tile boundaries, is only valid when the distance from each spiral core to the tile boundary is sufficiently large. Bohr *et al.*^{18,19} showed that for CGLE the tile boundaries are segments of hyperbolas with the two nearest spiral cores serving as foci and that the approximation is valid even when the separation between the cores is as small as one wavelength *λ*. The hyperbolic solution, however, only applies to spirals of opposite chirality (i.e., counter-rotating spiral waves).

A more general equation for the tile boundary between spirals of any chirality was derived by Zhan *et al.*^{27} Let the origins of the two spirals be **x** and $x\u2032$, their chiralities *σ*, $\sigma \u2032=\xb11,\u2009R=x\u2032\u2212x$, and let $x+r$ define a point on the boundary. Then the boundary is given by the solution to the differential equation

where $m=R/\lambda ,\u2009r\u2032=R2+r2\u22122rR\u2009cos\u2009\phi $, and $\phi $ is the angle between **r** and **R**. This equation can be solved numerically given an initial condition that lies on the line connecting the origins of the two spirals. Equation (10) was shown to accurately capture the tile boundaries not only for weakly nonlinear oscillations found in CGLE, but also for the Barkley model^{28} which, like the Karma model, describes an excitable medium supporting strongly nonlinear oscillations.^{29}

We checked the validity of this equation for the modified Karma model by comparing the analytic solutions to those computed using the cycle area method for the co- and counter-rotating phase-shifted solutions from Figs. 6 and 7, for which *m* = 3.18. Excellent agreement was found in all cases. Four examples with the analytic solutions superimposed on the cycle area plots are shown in Fig. 12. Numerical and analytical solutions also agree well for the three-spiral solution, as Fig. 9(c) illustrates.

It should be noted that both methods of computing the tile boundaries have advantages and drawbacks. Numerical solutions based on cycle areas do not require any assumptions (e.g., the Archimedean approximation for the phase) and can be computed in real time. On the other hand, they are only updated once per period (i.e., upon crossing of the Poincaré section), which may not be adequate for quickly drifting spirals. Furthermore, the shocks often do not entirely enclose each spiral (cf. Figs. 4(g)–4(h)); an additional construction is needed to form a closed boundary or determine the precise position of the boundary based on the transverse profile of the shock. Analytical solutions form closed boundaries, but require identification of the position of the cores and the initial condition (e.g., the point where the boundary crosses the straight line connecting the two cores). Analytical construction also requires an algorithm for determining the endpoints of each smooth segment of the boundary where three (or more) different tiles meet.

### B. Assembling a global solution

For computational domains with no-flux boundary conditions, shocks form naturally along the boundaries. Additionally, empirical observations for states with tile boundaries evolving slowly compared with the spiral rotation demonstrate that the level sets of both *v* (cf. Figs. 4(g)–4(h), 9(c), and 12) and *u* (not shown) are orthogonal to the shocks. Hence, the single-spiral solution on each of the tiles satisfies the no-flux boundary conditions on its entire boundary, whether it is external or internal with respect to the computational domain. We can therefore reduce the interaction between different spirals to the effect of boundary conditions, which only affects the solution in the interior of the tile through its shape (and the dynamics of the tile boundaries, if any). The effect of the boundaries on the enclosed spiral wave should be qualitatively the same regardless of the tile geometry. In particular, the temporal period of each spiral wave should depend on the size of the tile. For example, the period of unstable single-spiral wave solutions decreases from the asymptotic value *T*_{0} as the domain size *L* decreases for the Karma model on square domains.^{17} For domains of size comparable to the smallest tiles in Figs. 4(g)–4(h), the period is estimated to decrease by $O(10\u22122T0)$. Such differences in the periods of different spirals would lead to a relative residual of $O(10\u22121)$, which is consistent with the residuals reported in Figs. 4(c) and 4(d).

It is well known^{30,31} that if the frequencies *ω*_{1} and *ω*_{2} of two neighboring spirals differ, the boundary between them moves with velocity

where $k1$ and $k2$ are the wave vectors the two spirals would have on an unbounded domain at the location of the tile boundary. This is a consequence of the phase continuity at the tile boundaries: the frequencies of the two spirals become equal in a reference frame moving with velocity **c**. In particular, small differences in the frequencies (equivalently, periods) of two neighboring spirals lead to a slow motion of the boundary.

As we have shown elsewhere,^{17} in the Karma model the spiral cores are repelled from no-flux boundaries at distances of $O(2\u2113c$), where $\u2113c$ is the size of the spiral core. This appears to be a general result, as repulsive interaction was also found in other excitable systems.^{32,33} Hence, when the distance between a spiral core and the nearest tile boundary decreases to about $O(2\u2113c)$, the spiral core starts to drift away from the boundary. This is the mechanism that forces the core inside the small tile in Fig. 9 to drift relative to the cores in big tiles. Such core drift is also unavoidable for multi-spiral solutions with tiles of different sizes and will lead to the growth of big spirals at the expense of small ones, which can be considered as one of the mechanisms contributing to maintenance of spiral turbulence.

The tile-based decomposition suggests a natural approach to constructing global multi-spiral solutions from single-spiral segments satisfying local Euclidean symmetries. The results obtained for single-spiral solutions^{17} suggest that the solution on each tile would be defined either by a periodic solution or a relative periodic solution. If the tile boundaries do not move significantly during a typical period *T*_{0}, we can compute the solutions locally on each tile using the weighted Newton-Krylov method.^{17} The algorithm, however, will need to be generalized in such a way that updates of the initial condition at each step of Newton iteration preserve the continuity of both *u* and *v* fields (or equivalently the phase and amplitude) on the tile boundaries. This constraint, however, is not expected to present any conceptual difficulties.

## VII. CONCLUSIONS

To summarize, we have applied numerical methods originally developed for fluid turbulence to search for the exact coherent structures that may form a skeleton for the spatiotemporally chaotic dynamics produced by a prototypical monodomain model of cardiac tissue. We showed that these methods, designed to identify recurrent solutions in the presence of global symmetries, fail rather spectacularly for an excitable reaction-diffusion system whose dynamics is characterized by local, rather than global Euclidean symmetries. The origin of the failure was traced to the weak correlation between the dynamics of individual spirals which underpin spiral turbulence. As a result of this weak correlation, typical multi-spiral states display recurrent dynamics locally, but not globally. Locally the dynamics can be represented, to numerical accuracy, by periodic or relative periodic solutions, but globally neither periodic nor relative periodic solutions play a dynamically important role. Nonchaotic unstable solutions embedded in the chaotic set would have a more complicated nature and require development of novel computational approaches.

We propose one such approach based on the decomposition of the computational domain into sub-domains, or tiles, that each supports one spiral wave. Over short time scales (before individual spiral waves are destroyed by local instabilities), the dynamics for each near-recurrent multi-spiral state can be decomposed into the dynamics *of the tiles* (relative motion of tiles and changes in their shape associated with the differences in the spiral frequencies) and the dynamics of individual spiral waves *on the tiles* subject to no-flux boundary conditions at the tile boundaries. In particular, the dynamics of spiral waves on the tiles would be described by periodic or relative periodic solutions which correspond, respectively, to pinned and drifting cores.

It should be emphasized that the formalism based on decomposition into tiles is only expected to describe spiral turbulence during relatively quiescent intervals when no breakups or mergers of different spiral waves occur. Breakups and mergers are driven, respectively, by the alternans and meandering instabilities of individuals spirals^{17,20} and involve birth or annihilation of pairs of spiral cores with opposite chirality and the associated changes in the number of tiles. During such events, the cores move on the time scale comparable to the rotation period (or even faster) and it may not be possible to even define the tiles on the entire domain. These relatively active intervals should not be described using the proposed formalism and likely correspond to heteroclinic connections between different exact coherent structures. This can be rephrased in terms of the of skew-decomposition of the dynamics^{34–37} in the vicinity of relative solutions induced by local symmetries: the quiescent intervals correspond to the evolution mainly along the group manifold, while active intervals correspond to motion mainly transverse to the group manifold.

Future work will focus on a numerical implementation of the procedure we have described for constructing global solutions composed of locally symmetry-reduced solutions. Specifically, the method for computing single-spiral wave solutions on square domains should be validated for tiles with an arbitrary shape. Furthermore, it remains to be determined how individual solutions should be recombined such that the phases and amplitudes of neighboring spirals match. If successful, this framework may provide valuable new insight not only into spatiotemporally chaotic dynamics of cardiac tissue, but also other systems that exhibit local symmetries.

The definition of locality is, of course, relative. In the model considered here, local symmetries survive on domains, or tiles, whose dimensions significantly exceed the characteristic correlation length $\u2113c$ defined by the spatial extent of the adjoint eigenmodes for spiral wave solutions.^{38,39} For excitable systems, such as the FitzHugh-Nagumo, Barkley, and Karma models, these eigenmodes decay exponentially, reflecting the lack of any long-range correlations. Short-range correlations, however, may not describe all cardiac tissue models. For instance, bidomain models^{40} also include an additional Poisson equation for the extracellular potential, generating long-range correlations. Similarly, long-range correlations can arise as a result of stretch-activated feedback.^{41} Investigation of the relation between symmetries and the structure of exact coherent structures in bidomain models is of particular interest both because they provide a more realistic description of cardiac tissue, compared with the monodomain models, and because of the analogy with fluid dynamics where long-range coupling is due to the pressure field, which also satisfies a Poisson equation.

## ACKNOWLEDGMENTS

The authors would like to thank Vadim Biktashev 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.