Adjoint eigenfunctions of temporally recurrent single-spiral solutions in a simple model of atrial fibrillation

Spiral and scroll waves are robust solutions of excitable and oscillatory media in, respectively, two and three spatial dimensions. Stability of these solutions is described by the spectra of the evolution operators obtained by linearizing the governing equations. Although the governing equations describing homogeneous and isotropic media respect Euclidean symmetry (on the $ℝ2$ plane), the particular wave solutions do not. As a result, the respective evolution operators are generically non-self-adjoint, so that their right (conventional) eigenfunctions are not mutually orthogonal and do not coincide with the left (adjoint) eigenfunctions.

The spectra of the evolution operators have a useful dynamical interpretation: the eigenvalues describe the temporal evolution of different perturbations, with positive (negative) real part corresponding to unstable (stable) modes, and the eigenfunctions describe the spatial structure of the perturbation mode growing or decaying at the rate defined by the corresponding eigenvalue, while the adjoint eigenfunctions describe the sensitivity of different modes to perturbations in the initial conditions.

Although non-trivial wave solutions do not respect the underlying Euclidean symmetry of the evolution equation, the symmetry manifests itself in the emergence of marginal modes (with zero eigenvalues) in the spectrum. In two dimensions the spectrum contains three marginal modes which correspond to the three continuous Euclidean symmetries: translation in the two directions spanning the plane and in-plane rotation.¹ The corresponding modes are known as Goldstone modes in Quantum Field Theory and Pattern Formation. The corresponding adjoint eigenfunctions have been termed response functions.² For stable spiral waves the Goldstone modes represent the dominant degrees of freedom and, in the presence of weak interactions, their evolution is naturally described in terms of the response functions.

The earliest work illustrating the role of response functions in the context of reaction-diffusion systems is due to Keener³ who investigated the dynamics of scroll wave filaments. Scroll waves with a straight untwisted filament can be unstable even if the corresponding two-dimensional spiral wave solution is stable: the bend or twist of the filament leads to self-interaction which causes transverse motion of the filament that, for small curvature and torsion, can be described with the help of the response functions.

The weak interaction picture also motivated investigations which aimed to describe the drift of stable spirals in two dimensions. Biktashev and Holden⁴ showed that experimental and numerical results describing the dynamics of a spiral wave in the presence of resonant forcing, perturbations of parameters, or interaction with the boundary can be understood using an empirical model containing three coupled ordinary differential equations (ODE) for the position of the wave core and the phase of the wave. In a subsequent paper, the same authors⁵ showed that the ODEs can be derived with the help of the three response functions.

The response functions themselves, however, were not computed for any reaction-diffusion system until much later. Biktasheva et al.² computed them for the spiral waves in the complex Ginzburg–Landau equation (CGLE), which describes the generic dynamics of a broad class of spatially extended systems close to the onset of the oscillatory instability, and Henry and Hakim^6,7 computed response functions for scroll waves in the Barkley model of excitable media. In both instances the calculations were performed in a co-rotating reference frame which transforms the time-dependent solutions (relative equilibria) into steady states (absolute equilibria). This reduction based on the equivalence of the temporal evolution and rotation has been employed for computing the adjoint eigenfunctions in all subsequent studies.

With the help of the computed response functions it was possible to check that there is not just qualitative, but also quantitative agreement between numerical simulations and the ODE-based model for the resonant drift of spiral waves in CGLE⁸ in the Eckhaus-stable parameter regime. All three response functions were found to decay exponentially with the distance to the spiral core for the CGLE, in agreement with the analytical prediction. Exponential localization was later found even in the Eckhaus-unstable parameter regime.⁹ 

Response functions can be used to describe the interaction of spiral waves not only with boundaries, but also with other spiral waves. Indeed, exponentially decaying interaction for spiral wave solutions of CGLE was predicted previously using the amplitude/phase equation formalism.^10–12 The analytical results obtained for CGLE, however, cannot be extended to strongly nonlinear waves in excitable media, so the response function formalism becomes the only tractable means of predicting the evolution of spiral waves in response to internal or external perturbations.

Quite interestingly, the response functions were found to be exponentially localized also for the Barkley,⁷ FitzHugh–Nagumo,¹³ Oregonator,¹⁴ and Beeler–Reuter–Pumir¹⁵ models, suggesting that, as a rule, the spiral core acts as the organizing center for the wave, although there are rare counter-examples¹⁶ such as the Mornev model.¹⁷ 

The exponential localization of the response functions enables quantitative description of the dynamics of the spiral core as a singular forced object. As Ref. 18 puts it, “spiral waves look like essentially nonlocalized objects but behave as effectively localized particles.” As a result, despite the dissipative nature of excitable media, one finds a wave-particle duality that is more akin to that found in Hamiltonian systems. For example, Langham and Barkley^19,20 used the response function formalism to show that core of a resonantly driven spiral in a bounded domain moves almost like a classical particle, although reflections from the “walls” are characterized by a strongly nonlinear relation between the incident and reflected angle.

When the spiral or scroll wave is unstable, its dynamics cannot be described in terms of the marginal degrees of freedom (i.e., Goldstone modes and response functions). Instead, one also must consider the evolution of all the unstable modes. The only relevant study that we are aware of is due to Allexandre and Otani²¹ who considered the problem of feedback control of unstable spiral wave solutions of the Fenton–Karma²² model. In addition to the response functions, the eigenfunctions adjoint to all of the unstable modes were computed. The structure of the unstable adjoint eigenfunctions is especially important in the control problem as it allows significant optimization.^23,24

In all of the above studies only one type of spiral wave solution was considered—those described by relative equilibria—which can be reduced to a steady state in a rotating reference frame. However, most non-trivial applications involve more complex types of solutions. For instance, the simplest type of a spiral wave on a bounded domain of generic shape is described by a periodic solution. The simplest meandering spirals, even in $ℝ2$⁠, are described by relative periodic solutions that are reducible to periodic solutions in a translating and rotating reference frame. The generalization of even the most basic results obtained for relative equilibria to more complex types of solutions is far from straightforward. Furthermore, it is not understood what kind of effects are introduced by the intrinsic time-dependence of the shape of a spiral wave.

In this work we present a generic method for computing the Floquet spectra, including both the right and left eigenfunctions, of essentially time-dependent spiral wave solutions, such as absolute and relative periodic orbits, on domains of generic shape. Unlike methods based on symmetry reduction, the presented approach allows one to compute all the leading modes in the spectrum, not just the Goldstone modes and the corresponding response functions, for both stable and unstable solutions. Sec. II gives the mathematical preliminaries including the description of the mathematical model. Sec. III describes the properties of the Floquet eigenfunctions and their adjoints for drifting spiral waves described by (generalized) relative periodic solutions. The localization of adjoint eigenfunctions and the interaction of spiral waves with physical boundaries are discussed in Sec. IV. We end with some general conclusions and summary in Sec. V.

Our focus here is on the unstable spiral wave solutions that describe features of complex cardiac arrhythmias, such as fibrillation. Cardiac tissue is an example of excitable media, which are commonly modeled using reaction-diffusion equations

The interpretation of the dynamical variables $u=[u1…ul](t,x)$ depends on the details of the system under consideration. In monodomain models of cardiac tissue these represent the transmembrane voltage, ionic fluxes, etc. The diffusion tensor σ describes the intercellular couplings and the nonlinear function $f(·)$ describes the cellular dynamics. We use a modified form of the Karma model^25,26 with l = 2 and

where u₁ is the non-dimensional transmembrane potential and u₂ is a gating variable. This system of equations was formulated as a minimal model of cardiac excitation dynamics exhibiting alternans—the instability that is believed to be responsible for initiating and sustaining fibrillation.^27–29 Alternans manifests as a period-doubling bifurcation in isolated cells³⁰ and a Hopf bifurcation in tissue,³¹ where the temporal period and wavelength of the perturbation are nearly double those of the base excitation wave. The parameters of the model are as in Ref. 32, with the exception of the width of the switching function, which we set to s = 1.2571 throughout.

For $σ=D=const$⁠, Equation (1) with $x∈ℝ2$ and $t∈ℝ$ is equivariant with respect to continuous group actions from E(1) × E(2), consisting of temporal shifts $t→t′=t+τ$ and continuous translation and rotation transformations of the spatial coordinates $x→x′=Rϕx+h$⁠, where $Rϕ=exp(ϕ∂θ)$ is the operator of rotation by angle $ϕ$ and $∂θ$ is the generator of rotations. A consequence of this equivariance is that each solution lies on a manifold of equivalent solutions, generated by the action of symmetry transformations on the state.¹ In particular, the Goldstone Modes correspond to the generators of translational (⁠$∂xu$⁠) and rotational (⁠$∂θu$⁠) symmetries in space and forward-translation symmetry in time (⁠$∂tu$⁠) that span the tangent space of the group manifold.³³ 

The presence of these symmetries, in turn, allows the special types of traveling wave solutions (i.e., plane and spiral waves) described by relative equilibria. In particular, the rigidly rotating spiral waves which satisfy

can be found even on bounded domains as long as their shape is consistent with the rotational symmetry (e.g., circular domains).³⁴ The equivalence between rotation and time-translation explicit in (3) is the special property that enabled computation of the marginal modes by transforming the equation into a reference frame co-rotating with angular velocity ω, which transforms a relative equilibrium into an absolute equilibrium. This equivalence is also the reason there are three, rather than four, independent Goldstone modes for relative equilibria.

Although circular domains are computationally and mathematically convenient, they are physiologically irrelevant. In particular, cardiac tissue is not circular, and in more complex regimes multi-spiral states partition the domain into tiles, each of which supports a single spiral,³⁵ that never, even transiently, acquire a circular shape. Once the boundary loses circular symmetry, spiral wave solutions cannot be described by relative equilibria. Instead, one finds solutions such as pinned spirals (absolute periodic orbits)

or drifting spirals

where $UT$ is the time-T evolution operator, and $Th=exp(h·∇)$ is the translation operator. Relative periodic orbits on $ℝ2$ satisfy (5) exactly and become periodic solutions in a reference frame moving with velocity $c=h/T$⁠. On bounded domains, these solutions are elevated to generalized relative periodic orbits which are numerically exact (i.e., exact to the level of numerical accuracy) only in the interior of the domain, outside of a narrow boundary layer.³² On bounded domains, the time-dependence of spiral wave solutions is essential and cannot be removed by any symmetry transformation, requiring a different approach for computing the spectrum.

In this work we will assume that the spatial domain Ω is finite and has a square shape, and impose the typical for cardiac models “no-flux” boundary conditions

where $j=σ∇u$ is the flux and $n$ is the outward-oriented normal to the boundary $∂Ω$ of the domain. This choice represents a generic situation where the boundary conditions explicitly break the translational and rotational symmetry. The square shape is also representative of the tiling decomposition of multispiral solutions, where each tile is bounded by a small number of intersecting smooth curves meeting at finite angles.³⁵ 

The evolution Equation (1) along with boundary conditions (6) defines a solution

to the initial value problem, where $Ut$ is the time-t evolution operator of the nonlinear reaction-diffusion system (1). The temporal evolution of small perturbations $v(t,x)$ about this solution is described by the linearization of (1), with the boundary conditions $n·∇v(t,x)=0$ inherited from (6). The resulting (forward) evolution equation in the tangent space is

where $L[u]=D∇2+f′(u)$ is the instantaneous Jacobian operator and $f′(u)=∂f/∂u$ describes the variation of the cellular dynamics with respect to the state variables.

The propagator $Vt$ maps an initial infinitesimal perturbation $v(0,x)$ in the tangent space to a later time, $v(t,x)=Vtv(0,x)$⁠, where

Similarly, the action of the adjoint propagator $Vt†$ maps an infinitesimal perturbation $w(t,x)$ in the adjoint space to an earlier time, $w(0,x)=Vt†w(t,x)$⁠, where

and $L†[u(t,x)]$ is the adjoint of the instantaneous Jacobian operator (8). It is easy to check that $Vt†$ defines a time-reversed flow in the adjoint space

For T-periodic solutions, the eigenfunctions and eigenvalues of $VT$ satisfy

Correspondingly

where the asterisk denotes the complex conjugate. In other words $wi(T,x)$ and $vi(0,x)$ represent the left and right eigenfunctions of $VT$ corresponding to the eigenvalue Λ_i. Since $VT$ is generally non-self-adjoint, $wi(T,x)≠vi(0,x)$⁠, however the two sets of eigenfunctions are mutually orthogonal. The orthogonality relation can be written in a more useful and general form

for $0≤t≤T$⁠, where $vi(t,x)$ is defined as the solution of (8) that coincides with the right eigenfunction $vi(0,x)$ at t = 0 and $wj(t,x)$ as the solution of (11) that coincides with the left eigenfunction $wj(T,x)$ at t = T. For a time-periodic solution $u(t,x)$ with period T, $vi(t,x)$ and $wj(t,x)$ correspond to the Floquet modes.

Although the vast majority of studies of infinite-dimensional systems focus exclusively on the right (conventional) eigenfunctions, the importance of the left (adjoint) eigenfunctions is hard to overestimate. Indeed, an accurate estimate of the evolution operator $VT$ is crucial for a number of applications such as computing unstable solutions,³² feedback control,^21,23,24 and adjoint-based optimization. For infinite-dimensional systems, explicit evaluation of $VT$ is prohibitively expensive (or simply impossible). Such estimates are often constructed by using Arnoldi iterations³⁶ which generate approximations of the leading right eigenfunctions, but generally do not approximate the left eigenfunctions. An optimal estimate based on the truncated spectral decomposition

obviously requires both sets of eigenfunctions. In comparison, an Arnoldi-based approximation

e.g., used in GMRES methods,³⁷ uses synthetic adjoints $w¯i$ which satisfy (14) only in the M-dimensional Krylov subspace spanned by ${v1,…,vM}$ (i.e., $〈w¯j|vi〉≠δij$ if either i or j > M) and requires $M≫N$ to achieve a similar level of accuracy. Similarly, the adjoints are required to compute the coordinates a_i of disturbances in the tangent space

In a previous study³² we have introduced a numerical method which employs a Newton–Krylov solver for computing both the unstable generalized relative periodic orbits of (1) and their right eigenfunctions on bounded domains of arbitrary shape that break all the continuous symmetries of the underlying evolution equations. A snapshot of a single-spiral solution described by a generalized relative periodic orbit is presented in Fig. 1(a). This spiral wave solution has wavelength λ = 78 (2.04 cm) and period T = 54.74 (136.9 ms) on a domain of size L = 192 (5.03 cm). The unit conversion uses the typical size of a cardiac cell (262 μm) for nondimensionalization. When the tip $xo$ of the spiral (defined by $∂tu(0,xo)=0$⁠) is placed within 10 μm of the center of the domain, it drifts by just $|h|=O(10−11λ)$ over the course of the rotation, reflecting the discrete rotational symmetry of the problem associated with this particular initial condition. Hence, although it is a member of a class of relative periodic orbits, this particular solution is, to numerical precision, simply a periodic orbit. For other choices of x_o, we found that $|h|∝e−ζ/ℓc$⁠, where ζ is the distance from x_o to the nearest boundary and $ℓc$ is a numerically determined critical length scale.

This study extends the results of Ref. 32 by using a method (detailed in the  Appendix) that does not rely on a transformation to the co-moving frame, thus allowing computation of adjoint eigenfunctions for spiral wave solutions with arbitrary symmetry properties and temporal dependence. The spectrum of Floquet multipliers Λ_i for the spiral wave solution is shown in Fig. 1(b). The Floquet multipliers correspond to the eigenvalues of $VT$⁠; the corresponding Floquet exponents σ_i can be computed using the relation $Λi=eσiT$⁠. The spectrum is seen to include both a discrete and a continuous part. There are four unstable eigenvalues (two from the discrete and two from the continuous part) as well as a triply degenerate marginal eigenvalue (Λ = 1) associated with infinitesimal spatial and temporal translations, just like for relative equilibria.¹ The two Goldstone modes associated with infinitesimal spatial translations persist on bounded domains^1,32 due to the local translational invariance of (1) even though the global translational symmetry is broken by the boundary conditions (6).

We have computed the leading 130 left and right eigenmodes with high accuracy (see  Appendix). The movies of $vi(t,x)$ and $wj(t,x)$ for several dominant modes are included in the supplementary material (Sec. VI). As Fig. 1(b) shows, the eigenvalues associated with the left eigenfunctions are just as accurate as those associated with the right eigenfunctions. In particular, the eigenvalues associated with the marginal modes deviate from unity less than $O(10−7)$ for both the Goldstone modes and the response functions.

The Goldstone modes associated with translational symmetries are shown in Figs. 2(a) and 2(b) and 2(d) and 2(e). These correspond to the spatial derivatives of the initial condition $n·∇u$ along two orthogonal directions $n$⁠. Figures 2(c) and 2(f) show the Goldstone mode associated with temporal translation which corresponds to the temporal derivative $∂tu$⁠. Figures 2(g)–2(l) show the corresponding response functions. The orthogonality condition (14) does not completely fix the normalization of the two sets of eigenfunctions. As the eigenfunctions represent solutions to a linear problem (and are thus scale-independent), we are free to choose the absolute scale of each set, provided (14) is satisfied. Hence, we added an additional constraint $〈vi|vi〉=〈wi|wi〉$⁠, so that the scales of the dominant components of $vi$ and $wi$ are comparable.

As Fig. 2 illustrates, for the Karma model the second component (gating variable) of the right eigenfunctions is very small compared to the first component (voltage variable), $||v2i||=O(10−2||v1i||)$⁠, while for the left eigenfunctions the opposite is true, $||w1i||=O(10−2||w2i||)$⁠. This disparity is due to the difference in the time scales of cellular kinetics described by the nonlinear functions $f=(f̃1,ϵf̃2)$⁠, where $ϵ=0.01$ is a fixed parameter and both functions $f̃1$ and $f̃2$ and their partial derivatives are all O(1). For instance, if we rescale the two components of the right eigenfunction as $vi=(ṽ1i,ϵṽ2i)$⁠, then the evolution equation (8) can be rewritten as

where $∂f̃i/∂uj=O(1)$ for all i and j. Similarly, rescaling the components of the left eigenfunctions $w=(ϵw̃1i,w̃2i)$⁠, we can rewrite the evolution equation (11) as

Equations (18) and (19) have solutions both components of which are of the same order of magnitude, so we should indeed expect $||v2i||=O(ϵ)||v1i||$ and $||w1i||=O(ϵ)||w2i||$ for all eigenfunctions. Hence, in the remainder of the paper we will focus mostly on their dominant (unscaled) components $v1i$ and $w2i$⁠.

The most salient feature of the response functions in the Karma model is that they are very strongly localized near the core of the spiral, just like in most other excitable systems and the CGLE. To quantify this spatial localization we defined the amplitude³⁸ 

where r is the distance from $xo$⁠. The amplitude can be defined in a similar manner for all eigenfunctions, both left and right. As Fig. 3 shows, the response functions decay exponentially with r, while the amplitude of the Goldstone Modes, as expected, remains a constant outside of the core region. For comparison, the dashed line shows the exponential decay $A∝e−r/ℓc$ predicted numerically³² based on the scaling results for the spatial drift $h$ and the period T of the spiral wave, where $ℓc≈λ/20$ for the value of s used here. The spatial decay rate determined directly from the response functions gives a very close value $ℓc=0.0478λ≈λ/21$⁠, which confirms the conjecture³² that the scaling of $h$ and T is indeed controlled by the spatial structure of the response functions corresponding, respectively, to the spatial and temporal Goldstone modes.

The response functions do not decay with r monotonically. Instead, the dominant component w₂ displays pronounced oscillations, with roughly the same distance (⁠$Δr≈3.04ℓc$⁠) between the nodal lines for all three response functions. For the temporal response function the nodal lines form closed loops in the plane. In contrast, the nodal lines of the translational response functions form spirals (they are not closed curves). Correspondingly, the angular averaging destroys the underlying oscillation of the amplitude A(r) in the latter case, while in the former case the amplitude clearly shows the modulation superimposed on top of the exponential profile.

The method introduced here allowed us to compute not only the marginal eigenfunctions, but an entire spectrum of leading modes. In particular, two pairs of complex conjugate modes—the most unstable pair and a stable pair from the discrete part of the spectrum—are shown in Figure 4 and the corresponding angle-averaged amplitudes are plotted in Fig. 5. Again we find that the adjoint eigenfunctions are strongly localized in the core region. Their amplitude decreases exponentially, $A(r)∝e−r/ℓ−$⁠, just as it does for the response functions, albeit more slowly: the corresponding length scales are $ℓ−=1.72ℓc$ and $ℓ−=2.20ℓc$ for the unstable and stable modes, respectively. The right eigenfunctions show the opposite trend, their amplitude increases exponentially, $A(r)∝er/ℓ+$⁠. In particular, the stable mode is localized near the boundary, with amplitude growing on the length scale $ℓ+=7.04ℓc$⁠.

The length scale $ℓ+≈5λ$ for the unstable modes is extremely large, so they appear to be evenly distributed throughout the entire domain (cf. Figs. 4(a) and 4(b)). The shape of the modes clearly indicates that they describe the alternans instability which is characterized by the variation in the width of the excitation wave. This particular pair of modes describes discordant alternans: as Fig. 4(a) illustrates, at t = 0 the thickness of the excitation wave increases in some regions (positive values of $Re(v1)$⁠) and decreases in others (negative values of $Re(v1)$⁠), which corresponds to an increase (respectively, decrease) in the action potential duration (APD). The corresponding Floquet multipliers are complex (⁠$arg(Λ)≈±2π/3$⁠), rather than real and negative, as would be the case for a period-doubling bifurcation.

This is, however, not the only unstable mode. There are both absolutely unstable modes (characterized by $|Λ|>1$ or $Re(σ)>0$⁠) and convectively unstable modes (for which $Re(σ)<0$⁠) in the continuous part of the spectrum. A right eigenfunction with amplitude $A(r)>Cer/ℓ+$⁠, where C and $ℓ+$ are some positive constants, describes a convectively unstable mode, if $−c/ℓ+<Re(σ)<0$⁠, where $c=λ/T$ is the asymptotic conductive velocity of the spiral wave. Equivalently, the temporal decay $Re(σ)T$ must overwhelm the spatial growth $λ/ℓ+$ in time T and length λ for the mode to be convectively stable. Convectively unstable modes would produce a noticeable distortion of the spiral wave on sufficiently large domains. Unexpectedly, most of the absolutely stable modes computed for the solution shown in Fig. 1 grow with radial distance from the core so quickly that they are convectively unstable. As Fig. 6 illustrates, it is only a subset of the strongly contracting modes (⁠$|Λ|<O(10−1)$⁠) associated with relatively featureless eigenfunctions (i.e., for which $ℓ+=O(λ)$⁠) which are both convectively and absolutely stable. Convective instabilities and exponential growth of eigenfunctions for defect-modulated waves in reaction-diffusion systems have been investigated previously in one³⁹ and two^40,41 spatial dimensions, and their role is reasonably well-understood.

As Fig. 1 shows, the continuous spectrum crosses the unit circle near $Λ=e±5πi/6$⁠, indicating that on $ℝ2$⁠, one would expect to find an infinite number of modes close to the boundary of absolute instability. For the relatively small domain size considered here (⁠$L=2.46λ$⁠), only a single complex pair of modes from the continuous spectrum is absolutely unstable (cf. Figs. 7(a)–7(d)). However, the continuous spectrum contains a large number of modes that are convectively unstable (two examples are shown in Figs. 7(e)–7(l)). Figures 7 and 8 show the eigenfunctions of three leading modes from the continuous spectrum and their amplitude.

The leading modes in the continuous spectrum exhibit spatial localization trends similar to those from the discrete spectrum. In particular, the left eigenfunctions are localized in the core region, while the right eigenfunctions are localized near the boundaries (cf. Fig. 7). The mode amplitudes, however, are not given by pure exponentials, but rather a product of an exponential and a power, i.e.,

with $α≈2$ and $ℓ−≈2.5ℓc$ for the left modes, and $ℓ+≈16ℓc$ for the right modes, with $Λ±=−0.88±0.48i,−0.83±0.50i$⁠, and $−0.81±0.50i$⁠. The spatial structure of the right eigenfunctions suggests that these modes are also related to alternans, although the width variation is mixed with bending of the excitation wave (in some regions both the leading and the trailing shock are displaced in the same direction, rather than in the opposite directions, as would be the case for alternans).

To sum up, for the domain size considered here, there are several alternans modes in the Karma model. Classical alternans (cf. Figs. 4(a) and 4(b)), the dominant mode of instability, lacks spatial localization, while there is another unstable alternans mode that is localized near the boundary (cf. Figs. 7(a) and 7(b)). The adjoints of these modes are all strongly localized near the core. In comparison, in the three-variable Fenton–Karma model,²² the alternans modes of an unstable spiral wave computed on a disk of comparable size²¹ are characterized by adjoint eigenfunctions that show almost no attenuation with r. Our results show that, unlike the Fenton–Karma model where the development of alternans appears to be sensitive to perturbations over the entire domain, in the Karma model the development of alternans is only sensitive to perturbations near the spiral core.

In the conclusion of this section, we discuss the structure of some of the (absolutely) stable modes from the continuous part of the spectrum. The dominant components of the modes with Floquet multipliers $0.4<|Λ|<0.6$ are shown in Fig. 9. Similar to the unstable and weakly stable modes, the left eigenfunctions are found to be localized near the core of the spiral and the right eigenfunctions near the domain boundary, although the localization is weaker than for the modes with larger $|Λ|$⁠. The amplitudes (cf. Fig. 10) are again found to increase/decrease exponentially (aside from some weak modulation) with r on length scales $ℓ±=O(6ℓc)$⁠.

Although spatial localization of adjoint eigenfunctions appears to be an almost universal property of spiral wave solutions, it has only been understood for the response functions associated with a spiral wave solution in the CGLE.^2,9 Extending these results for generic excitable systems has proved difficult due to the strong nonlinearity of the evolution equations. However, we can make progress in certain limits. Although the spiral wave solution investigated here formally corresponds to a relative periodic orbit, it is nearly indistinguishable from a rigidly rotating spiral wave (i.e., relative equilibrium) inside a circle of radius $L/2$⁠. On an infinite domain, far from the origin, the Archimedian approximation applies

where $ξ=r+lθ−ct, l=λ/(2π), c=λ/T$⁠, and $u0(ξ)$ is periodic with period λ. The adjoints can then be written in the form $w̃(r,θ,t)=w¯(ξ)eimθeγt$ and the second of the two equations in (19) becomes

For $r≫λ$⁠, the curvature of the spiral wave can be ignored, and (23) simplifies, yielding

Since the partial derivatives $∂f̃i/∂uj$ depend only on $u0$⁠, they are periodic in ξ and, according to the Floquet theorem, Equation (24) has solutions in the form $w¯2i(ξ)=ŵ2i(ξ)ekiξ$⁠, where $ŵ2i(ξ+λ)=ŵ2i(ξ)$⁠. We therefore should expect the adjoints to grow or decay exponentially with ξ or r as long as $Re(ki)≠0$⁠. This result generalizes the prediction of exponential far-field dependence of left eigenfunctions in one spatial dimension.^40,42 Note that $γ=σi*+cki$⁠, where

is the corresponding Floquet exponent.

We can make further progress in various special cases. For instance, when $ϵ≪|σi|$ (e.g., for strongly stable modes), we have $ŵ2i=const$ and $σi*=−D22ki2$ with solutions $ki=|ki|eiχi$⁠, where

and n_i is an integer. Therefore solutions $w2i$ can both decay and grow exponentially with r on finite domains. On an infinite domain, however, only solutions that decay exponentially with r are allowed (⁠$Re(ki)<0$⁠), which explains exponential localization of the adjoints with a length scale $ℓ−=Re(ki)−1∼D22/|σi|$⁠.

In the limit $D22→0$ (which is the typical case considered in models of cardiac tissue), D₁₁ becomes the only parameter in Equation (11) with the dimension of length. Hence the localization length scale for all slow (unstable, marginal, and weakly stable) modes characterized by the time scale $ω−1$ can be found using dimensional analysis, which yields $ℓ−∼D11/ω=5.90$⁠. This is fairly close to the numerical value found for the response functions, $ℓc≈3.72$⁠, for $D22/D11=0.05$⁠. Indeed, this is not entirely unexpected: as Fig. 11 illustrates, the localization length scales for the amplitude of all three marginal adjoint eigenfunctions (⁠$wx, wy, wt$⁠) depends rather weakly on D₂₂.

For strongly stable modes the relevant time scale can be quite different, although this difference may only become apparent for very quickly decaying modes that are not resolved in the numerics. For instance, for $|Λ|=0.01$ we have $|σi|∼−ln|Λi|/T=0.08$ and therefore $ℓ−∼D22/|σi|=1.6$⁠, which is of the same order of magnitude as the value we found for the marginal modes. The similarity of the length scales predicted for slow modes and the strongly contracting modes explains the relatively small variation in the localization between adjoint eigenfunctions throughout most of the spectrum.

It is also instructive to compare the structure of the response functions adjoint to the Goldstone modes associated with spatial translations with the structure of the shift map that defines how interaction with a no-flux boundary affects the drift of relative periodic solutions. We have shown previously³² that the distance ζ_n between the spiral core and the (planar) boundary after n periods of the revolution can be described by a map

The roots $0<ζ0<ζ1<⋯$ of the shift function $hn(ζ)$ define the equilibrium separation values. When the distance between the origin $xo$ of the spiral and the closest boundary is equal to one of these equilibrium values, the spiral wave will drift tangentially to the boundary. An equilibrium ζ_k is stable provided $|1+h′n(ζk)|<1$ and unstable otherwise. Since $|h′n(ζ)|≪1$⁠, this inequality is equivalent to $h′n(ζk)<0$⁠. In particular, $ζ0=5.36ℓc$ is a stable fixed point, while $ζ1=7.85ℓc$ is unstable. The existence of stable equilibria suggests the presence of bound states, where a spiral would drift along a planar no-flux boundary forever. Similar bound states were found for resonantly driven spirals next to effective boundaries.^19,20

The drift of spiral waves caused by interaction with the boundaries can be understood by considering the relation between solutions of (1) on bounded domain and on $ℝ2$⁠. Consider the flux $j=σ∇u$⁠, where the diffusion tensor $σ=D$ inside a bounded domain Ω and σ = 0 outside. In this case, the no-flux boundary condition (6) is equivalent to the inclusion of an additional term

on the right-hand-side of (1) defined on $ℝ2$⁠. Here $δ∂Ω$ denotes a one-dimensional delta function localized at $∂Ω$⁠, such that in a small neighborhood of every point $xb∈∂Ω$

In the absence of this additional term, (1) possesses a spiral wave solution rigidly rotating around the tip x_o, which corresponds to the relative equilibrium on $ℝ2$⁠. The introduction of this term generates a perturbation to the dynamics of all the modes of this solution. In particular, the perturbation along the Goldstone modes $vx, vy$⁠, and $vt$ will generate, respectively, the drift of the spiral core in the $x̂$ and $ŷ$ directions and a phase shift (rotation). It will be convenient for us to define $vq(t,x)$ satisfying (8) such that

where $q={x,y,t}$ and

where the coefficients α_qj are chosen such that the orthogonality condition $〈wp|vq〉=δpq$ is satisfied at t = T. Thus defined, $wq(t,x)$ will satisfy (11) and be orthogonal with respect to $vq(t,x)$ at all t. The spatial drift is then given by⁴³ 

It should be stressed that this relation is only exact for infinitesimal perturbations, while (28) is not infinitesimal. Nonetheless, (32) provides a fairly accurate description of the drift, as we will see below.

We can make further progress assuming the boundary $∂Ω$ is a smooth curve. With the help of (28) the components of (32) can be rewritten as

where $dl$ is the arclength element along the boundary. Placing the origin of the coordinate system at the tip of the spiral wave, we can write $wq(t,x)=ŵq(t,x)e−r/ℓc$⁠, where $r=x2+y2$ and $ŵq(t,x)=O(1)$⁠. Since $(n·∇)u(t,x)=O(1)$ as well, the integral is dominated by the region of the contour around the point x_b closest point to origin. Let us orient the coordinate axes such that $xb=(ζ,0)$ and $n=x̂$⁠. Since on the contour of integration

the integral can be evaluated using the saddle point method yielding

Higher order corrections can be easily generated and scale as $(ζℓc/λ2)m$ relative to (35) with integer $m≥1$⁠, where $λ2≫ζℓc$⁠, for this spiral wave solution. Outside of the core, the Archimedian approximation (22) can be used to evaluate the spatial derivative in (35)

where $ξb=ζ−ct$⁠, if we choose θ = 0 on the x axis.

The displacement $h=(hx,hy)$ can be found by integrating (35), where x_b will be a function of x_o. For a boundary with low curvature (⁠$κ≪ζ−1$⁠), we only need to keep track of the change in the distance ζ which satisfies

In particular, when $ζ≫ℓc$ we can neglect the change in ζ, so that the normal component of the shift $hn=n·h=hx$ is given by

The dependence of the drift function on the parameters ϵ, D₁₁, and D₂₂ can be made more explicit by factoring out the dependence of the components of the solution and the response function on the small parameter ϵ: $ŵx=(ϵw̃1x,w̃2x)$ and $u′0=(ũ′0,1,ϵũ′0,2)$⁠. Then (38) can be rewritten as

where

is a non-dimensional function

and $ℓr$ are characteristic length scales that represents the dependence of the amplitude of the response function $wx$ on parameters.

Rather predictably, we find that the dominant contribution to the drift is determined by the component of the state which has the largest diffusion constant. In the Karma model considered here $D11≫D22$⁠, and so it is the first component (which corresponds to the voltage variable) that controls the drift. For the inner product generally, where the contributions from the first and second components are unweighted, it is the contribution of the second component which dominates. Hence the normalization condition for $wx$ dictates that $w̃x∝ℓf/ℓc2$⁠, where $ℓf$ is the length scale on which the second component of $u0$ varies. Using the Karma model (1) it is straightforward to show that $ℓf=λ/(2ϵ)$ (cf. the spatial Goldstone modes in Figs. 2(d) and 2(e) and the second component of the solution shown in Fig. 1(a)). Therefore, we can define $ℓr=ℓc2/ℓf=2ϵℓc2/λ$ and

This length scale is quite different from a naïve guess $ℓd∼ℓc$ based solely on the scaling of the response functions.

As Fig. 12(b) shows, both components of $w̃x$ exhibit pronounced oscillatory dependence on the distance r from the origin. The scaled shift function $h¯n(ζ)$ inherits this oscillatory dependence: the expression (40) is shown as the dashed line in Fig. 12(a). The amplitude of the oscillation is nearly constant for $ζ/ℓc≫1$⁠, and to a high accuracy we can fit $h¯n(ζ)=A sin(πζ/Δζ)$ with constant A and $Δζ$⁠. The corresponding scaling for the shift function

yields a prediction which is in good agreement with both our previous numerical results³² and with the drift equation (33) integrated over one temporal period (cf. Fig. 12(a)). The fitting parameters for all three cases are given in Table I. In particular, we find that the values of A are O(1), which supports our choice of the length scale $ℓr$⁠. The values of $ℓc$ found by fitting the shifts correspond reasonably closely to the localization length scale of the response function. Moreover, the spacing $Δζ$ between the roots of the shift function corresponds well to the distance $Δr$ between the nodal lines of the response function. This confirms our conjecture³² that the interaction of spiral waves with a physical no-flux boundary is controlled by the response functions, just as in the case of resonantly driven spirals interacting with effective boundaries formed by a step-wise change in the excitability of the medium.^19,20

The saddle point approximation noticeably overestimates the magnitude of the shift due to interaction of the spiral wave with the boundary, while integrating (33) directly produces an estimate that is in reasonable quantitative agreement with the result of direct numerical simulations. While the saddle point approximation can be easily improved at larger $ζ/ℓc$ by retaining higher order terms in ζ, its main value is in uncovering the explicit dependence of the drift on various parameters of the problem. In practice, it is only the first root of the drift map, ζ₀, which is likely to play any role in the dynamics of spiral waves. The rest of the equilibria are essentially marginally stable, and the shift becomes exponentially small. The saddle point approximation produces a reasonably good prediction for the value of ζ₀. Even the prediction of the spacing $Δζ$ between the equilibria is in fairly good agreement: the error is only about $0.0485ζ0$⁠.

Qualitatively similar results were obtained for spiral interaction in CGLE using the amplitude equation formalism⁴⁴ and for spiral interaction with domain boundaries and defects in excitable systems using the kinematic theory.⁴⁵ In the latter case, however, the interaction strength was predicted to decay super-exponentially fast (i.e., as $exp(−ζ3)$⁠), a scaling our results do not support.

On domains of arbitrary shape, pinned and drifting spiral wave solutions of excitable systems are described, respectively, by temporally periodic and generalized relative periodic solutions. We have developed a general numerical procedure that allows computation of the leading adjoint eigenfunctions for unstable spiral wave solutions of such types. In particular, we computed hundreds of the dominant adjoint eigenfunctions for (slowly) drifting single-spiral wave solutions of the Karma model.

Just like for spiral wave solutions described by relative equilibria on circular domains, we found that the response functions, or adjoint eigenfunctions that correspond to marginal degrees of freedom, are exponentially localized in the vicinity of the spiral tip. The localization length scale of the response functions found numerically is in good agreement with the order-of-magnitude estimate $ℓc∼D11/ω$ based on dimensional analysis, where D₁₁ is the diffusion constant associated with the fast variable and ω is the angular frequency of the underlying spiral wave solution.

Adjoint eigenfunctions associated with other leading modes, both stable and unstable ones, were also found to be exponentially localized in the vicinity of the spiral tip, with the corresponding localization length scale $ℓ−$ larger than $ℓc$⁠. For strongly stable modes it can be shown more rigorously that $ℓ−∼D22/|σ|$⁠, where D₂₂ is the diffusion constant associated with the slow variable and σ is the corresponding Floquet exponent.

The significance of response functions for the dynamics of isolated spiral waves on $ℝ2$ is well understood.^9,13,38 The spatial and temporal response functions determine the effect of small perturbations in the initial conditions or the evolution equations on, respectively, the drift of the spiral and its rotation speed.^14,43 In particular, the response functions have been used to describe the interaction of spiral waves with tissue heterogeneties.^15,46,47 Our results further show that the spatial response functions also determine the interaction of spiral waves with physical no-flux boundaries and, by extension, the interaction with neighboring spirals through tile boundaries with effective no-flux boundary conditions.^32,35 Specifically, the spatial response functions $wx$ and $wy$ define the shift function $h(ζ)$ which describes the displacement of the spiral wave origin due to the interaction with the boundary over one temporal period.

The rest of the adjoint eigenfunctions have received very little attention in the literature dealing with excitable systems in general and the dynamics of cardiac tissue in particular. The spatiotemporal structure of unstable and weakly stable adjoints, however, is critically important for understanding spiral wave breakup and chaotic dynamics featuring multiple interacting unstable spiral waves. For instance, it is well known that stable spirals whose cores are sufficiently well separated can be considered effectively independent. The same is also true of unstable spirals.³² The spirals begin to interact at smaller separations, with interaction that can be conveniently described with the help of the adjoint eigenfunctions. Our results show that not only the position of the core and the phase of a spiral wave, but also its stability should be affected by neighboring spirals.

Furthermore, adjoints associated with slow modes play a crucial role in the design of feedback control methods aimed at suppressing spiral wave instabilities (such as alternans).²¹ The spatial localization of the adjoints associated with unstable modes indicates that feedback is most effective when it is applied close to the core of a spiral wave. Furthermore, the spatial alternation of the phase of the adjoints associated with all slow modes (unstable, marginal, and weakly stable) in the core region significantly attenuates the effectiveness of spatially uniform perturbations^4,19,48 on the dynamics of spiral waves. This suggests that spatially localized perturbation, e.g., those due to virtual electrodes^49,50 should be much more effective for control of spatiotemporally chaotic regimes, such as fibrillation.

See supplementary material for the temporal evolution of the leading eigenfunction pairs, as well as the rescaled response function used in the computation of the drift.

This material is based upon the 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.

The right eigenfunctions of the time-evolution operator $VT$ have been computed via Arnoldi iteration³⁶ which involves time-integration of the linearized equation (8). It should be pointed out that, since the instantaneous Jacobian L is a function of the reference state $u(t)$⁠, both (1) and (8) are time-integrated simultaneously to avoid storing and retrieving the reference solution. The same spatial (2nd order finite difference on a square mesh) and temporal (4th order fully explicit Runge–Kutta) discretization scheme is used for both equations.³² 

However, the same approach cannot be used to compute the left eigenfunctions, since (11) should be integrated backwards in time and (1) cannot be time-integrated in the reverse direction. Furthermore, as the evolution equation is quite stiff, fairly small time steps have to be used (⁠$O(104)$ time steps per period). One period of a fully resolved solution corresponds to about 8 GB of data, which may not fit in RAM. Hence, the entire reference solution $u(t)$ must be pre-computed, stored, and then retrieved during the time-integration of (11).

The spatial discretization and the time-stepping of the adjoint tangent evolution (11) are the same as those for (1) and (8). This choice was made out of necessity: the discrete adjoint of an explicit Runge–Kutta method is at least semi-implicit.^51,52 For a partial-differential equation, this requires solution of an infeasibly large linear system at every time step. Furthermore, the discretization of (11) is sufficiently precise that the solution of the large linear system is unnecessary.

Runge–Kutta integrators for (11) require evaluation of $u(t)$ at intermediate points between the time steps, while the solution $un=u(tn)$ is only known at discrete times $tn=nΔt$⁠. To preserve the accuracy of time-integration of (11), we use a high-order interpolation of $u(t)$⁠. That is, for an integrator of order $O(Δtp)$⁠, we use an interpolant uniformly accurate on the interval $t∈[tn,tn+1]$ to order $O(Δtq)$⁠, with $q≥p$⁠. Following the methodology of Enright et al.,⁵³ the $O(Δt4)$ interpolant for the classical Runge–Kutta method used in this work is

where $u3(t+ηΔt)$ is the value of the $O(Δt3)$ interpolated state at time $tn+ηΔt$⁠, which for this method has $η=1/3$⁠. The coefficients for the fourth order interpolant are

where $d4,1(τ)=1−d4,0(τ)$⁠. Due to interpolation, it takes roughly twice as long to integrate the adjoint evolution equation (11) compared with the tangent evolution equation (8) (e.g., 31 h vs. 70 h to compute the spectra shown in Fig. 13 on a single NVIDIA Tesla K20 GPU).

The Arnoldi iteration procedure we used to compute the spectrum (eigenfunctions and eigenvalues) of the finite-time tangent evolution map $VT$ can be applied to any time-periodic solution $u(t)$⁠. It relies on time-integration of (8) forward in time (with the final state $vi(T)$ for one iteration used to define the initial condition $vi+1(0)$ for the next iteration), to generate a sufficiently large orthonormal basis ${v̂1,…,v̂k}$ which defines the Krylov subspace that contains the dominant modes and an approximate Jacobian H_k in that basis

The adjoint spectrum (eigenfunctions and eigenvalues of $VT†$⁠) is computed using a similar approach (time-integrating (11) backwards in time with the final state $wi(0)$ for one iteration used to define the initial condition $wi+1(T)$ for the next iteration), yielding an orthonormal basis ${ŵ1,…,ŵk}$ and an approximate Jacobian $H′k$

The eigenvalues and eigenvectors of the matrices H_k, $H′k$ approximate the eigenvalues and (projected) eigenfunctions of the time-evolution operator and its adjoint, respectively.

The generation of eigenfunctions via the Arnoldi iteration relies on the closure of the reference solution $u(t)$⁠, as iteration forward in time by $VT$ requires both the beginning and endpoints of the trajectory to coincide in state-space. If the beginning and end states are distinct, then so are the associated tangent spaces, and successive applications of the map $VT$ become ill-defined. The set of basis vectors that defines the Krylov subspace is constructed by performing Gram–Schmidt orthogonalization of the sequence of Krylov vectors a_k—collocated in state space about the initial condition $u(0)$—constructed using the iterative relation

which spans the dominant subspace of the operator. For relative periodic orbits $u(T)=UTu(0)=g−1u(0)$⁠, where $g≠1$ is a symmetry transformation. This allows generalization of the iterative relation