Determining the source of period-doubling instabilities in spiral waves
DDetermining the source of period-doubling instabilities in spiralwaves
Stephanie Dodson ∗ Bj¨orn Sandstede † May 30, 2019
Abstract
Spiral wave patterns observed in models of cardiac arrhythmias and chemical oscillations de-velop alternans and stationary line defects, which can both be thought of as period-doublinginstabilities. These instabilities are observed on bounded domains, and may be caused by thespiral core, far-field asymptotics, or boundary conditions. Here, we introduce a methodology todisentangle the impacts of each region on the instabilities by analyzing spectral properties ofspiral waves and boundary sinks on bounded domains with appropriate boundary conditions.We apply our techniques to spirals formed in reaction-diffusion systems to investigate how andwhy alternans and line defects develop. Our results indicate that the mechanisms driving theseinstabilities are quite different; alternans are driven by the spiral core, whereas line defects ap-pear from boundary effects. Moreover, we find that the shape of the alternans eigenfunction isdue to the interaction of a point eigenvalue with curves of continuous spectra.
Keywords: spiral waves, reaction-diffusion systems, stability, alternans, period-doubling
AMS subject classifications:
Systems of oscillatory and excitable media frequently express spiral wave patterns. Spiral waves areobserved in laboratory settings in chemical oscillations in the Belousov-Zhabotinsky reaction [44,46]and cell signaling in slime molds [30], and have been associated with arrhythmic heart rhythms[28, 42, 43]. These systems support rigidly rotating spirals with constant shape, but transitions tocomplex dynamics and unstable spirals are common.In cardiac dynamics, accelerated tachycardiac rhythms have been linked to electrical activity orga-nized as rotating spiral wave patterns on the surface of the heart. The transition from tachycardiac ∗ Division of Applied Mathematics, Brown University, Providence, RI: stephanie [email protected] † Division of Applied Mathematics, Brown University, Providence, RI: bjorn [email protected] a r X i v : . [ m a t h . D S ] M a y o fibrillation is believed to be initiated by spiral wave breakup [27, 31]. Clinical studies indicate aprimary driver of breakup is conduction block following a long-short temporal modulation of theaction potential duration, in what is known as the alternans instability. Alternans are visible onelectrocardiograms and have become a clinical warning sign of sudden cardiac death [27, 31]. Inspiral waves, alternans physically corresponds to variation in spiral band width (Figure 1). For adetailed review of spiral waves in cardiac dynamics, see the review article [2] and recent resultspublished in the special issue [8]. Time (a) (b)
Figure 1: (a) Stationary line defect in the w -component of the R¨ossler system. (b) Timeevolution of alternans instability in the u -component of the Karma model on a square ofside length 16cm with homogeneous Neumann boundary conditions. System parameters asdefined in Section 2.Spirals are also produced and studied in chemical oscillations, for example the Belousov-Zhabotinskyreaction [44, 46]. In these systems, spirals are experimentally observed to form stationary line de-fects [26, 45] (Figure 1a), which have been reproduced in numerical simulations [16, 37]. Across thedefect lines, wave amplitudes are out of phase (Figure 1). Both alternans and line defects lead toa spiral wave with twice the period of the original planar wave.Nonlinear reaction-diffusion systems qualitatively capture transitions to complex meandering, drift-ing, and the period-doubled line defects and alternans patterns. In these systems, planar spiralwaves are stationary solutions in a rotating polar coordinate frame and converge to one-dimensionalperiodic travelling waves away from the core. Stability and bifurcations can be studied by con-sidering the spectra of the operator obtained by linearizing the nonlinear system about the spiralwave solution. The spectrum consists of isolated eigenvalues and a set determined by the operatorin the far-field limit.Bounded domains are of interest in applications to cardiac dynamics and laboratory experiments.Neumann boundary conditions naturally represent lower conductance tissue separating regions ofthe heart or the physical walls of containers. Mathematically, on finite domains planar spiralwaves are truncated and matched with a boundary sink, which adds extra structure to the spiralwave and alters the spectrum of the linear operator [13, 32, 33]. The boundary sink itself directlycontributes an additional set of eigenvalues, and the finite domain modifies the spectrum associatedwith the far-field dynamics. Furthermore, radial growth in the eigenfunctions is permitted, andthose that would not be integrable on the full plane now emerge as true eigenfunctions on the2ounded domain [13, 32, 33]. These eigenvalues are associated with intrinsic properties of thespiral wave and are attributed to the spiral core. The spectrum of the operator on a boundeddomain is therefore a union of three disjoint sets that are associated with instabilities from the far-field, boundary conditions, and core. Knowing which set unstable eigenvalues belong to providesinformation about how instabilities will manifest themselves on unbounded or bounded domains.Meander and drift instabilities are the result of a Hopf bifurcation originating from the core: theemerging dynamics is understood through actions of the symmetry group of translations and rota-tions on the plane and a center manifold reduction [5, 39]. However, previous studies investigatingalternans and line defects provide inconsistent and incomplete evidence for which spectral set theunstable eigenvalues belong to.Due to the clinical significance, the alternans instability has been a recent area of focus in thecardiac dynamics community. In single cells, alternans are widely attributed to a period-doublinginstability observed in simple 1D maps [18]. However, this condition, known as the restitutionhypothesis, has received contradictory evidence [7, 10] and does not appear to be relevant forexcitable tissues that support traveling waves. The formation and stability of alternans in wavespropagating on a ring and line have been analyzed with kinematic descriptions [12] and throughlinear stability analyses [3, 4, 6, 11]. Stability analysis in 1D predicts that alternans are the resultof a Hopf bifurcation [17], yet analysis of the 1D traveling waves cannot fully capture 2D features.Linear stability analysis of spirals on bounded domains in the Karma and Fenton-Karma modelsfound a variety of unstable eigenmodes responsible for the formation of alternans [1, 23, 24]. In [23]and [24], Marcotte and Grigoriev find that formation of alternans depends on the domain size.Furthermore, they determine that the alternans eigenmodes are not spatially localized near thecore.The rigorous analysis of spiral waves in [37] indicates period-doublings are initiated by a seriesof Hopf bifurcations with imaginary parts of the eigenvalues sitting robustly at multiples of halfthe spiral frequency. These Hopf eigenvalues may be induced by period-doubling of the far-fielddynamics or boundary sinks. Stationary line defects are hypothesized to stem from bifurcations ofthe boundary sink, but no direct evidence supporting this claim was found in [37] .The goal of this paper is to further investigate how and why these period-doubling like instabilitiesarise on bounded domains. Specifically, we seek to answer which spectral set the unstable eigenval-ues belong to and gain a better understanding of how the spirals destabilize and on what domainsthe instabilities are relevant. To tackle this problem, we introduce a methodology for disentanglingthe contributions of each region by forming related patterns on domains whose spectra will containeigenvalues arising from a subset of the resulting spectra. Three cases are considered and comparedwith essential and absolute spectra from wave trains: (1) a spiral on a bounded disk with Neumannboundary conditions to provide the full spectrum, (2) a boundary sink to demonstrate effects ofboundary conditions, and (3) a spiral on a disk with non-reflecting boundary conditions to removeboundary eigenvalues. 3e apply this methodology to reaction-diffusion models and discover that the mechanisms drivingalternans and line defects are rather different. We find that line defects arise from the boundarysink and thus will only appear under the correct conditions on bounded domains. In contrast,alternans originate from instabilities associated with the spiral core and will develop independentof the domain. Furthermore, the spectral computations reveal that the structure of the alternansinstability develops due to an interaction of an unstable point eigenvalue and curves of continuousspectra. Our results have important consequences for reproducing patterns such as line defects andprovide justification for extending analysis of alternans from simple bounded disks to the complexgeometry of the heart.In the sections that follow, we begin by describing the mathematical set-up and notation of thereaction-diffusion models. We include a review of relevant spectral properties for operators on theplane and how these properties are modified by bounded domains. Procedures used to computethe spirals and spectra are described in the methods section. Finally, we present the results of theanalysis applied to the R¨ossler and Karma models to study line defects and alternans, respectively. Reaction-diffusion systems display a rich set of patterns and are commonly used to model systemsin biology and nature. General planar reaction-diffusion systems are of the form U t = D ∆ U + F ( U ) , U ∈ R n , D ∈ R n × n , x ∈ R , (1)where U = ( u , . . . , u n ) T is a vector of species that diffuse at rates given by the nonnegative elements δ i of the diagonal matrix D , and ∆ is the Laplace operator. Kinetic reactions of the different speciesare captured by the typically nonlinear function F ( U ).Cardiac models range in complexity from biophysically detailed ion-channel models to simplifiedsystems which capture qualitative features, with both categories falling under the reaction-diffusionframework. The Karma system is a two-variable reduction of the Noble ion-channel model [25] andwas developed to be a simplified model that reproduces alternans [20, 21]. The model is given by u t = 1 . u + 400 (cid:18) − u + (cid:0) . − v (cid:1) (1 − tanh( u − u (cid:19) (2) v t = 0 . v + 4 (cid:18) − e − µ K ϑ ( u − − v (cid:19) , where the fast variable u represents voltage and v acts as a slower gating variable. As in [1, 23, 24],we use the function ϑ ( u ) = (1 + tanh(4 u )) /
2. Alternans are observed in this system when the realbifurcation parameter µ K is increased above one [21].The R¨ossler model is commonly used to study chaotic turbulence in chemical oscillations, and4s known to produce spirals with line defects. This three-variable system is also of the generalreaction-diffusion form and is given by u t = 0 . u − v − w (3) v t = 0 . v + u + 0 . vw t = 0 . w + uw − µ R w + 0 . . Bifurcations to line defects are observed for parameter values µ R > µ K and µ R , respectively. We remark that in (2) and (3) we selected values for several parameters that areoften allowed to vary. We refer to Table 1 in the appendix for the general form of these models. In our context, periodic traveling waves, also referred to as wave trains, serve as building blocks ofspiral waves. Therefore, we first consider the existence and stability properties of wave trains on R and their restriction to bounded domains before describing spiral waves on the plane and boundeddisks. Further details can be found in [13, 19, 32, 33]. On R , the reaction-diffusion system (1) reduces to U t = DU xx + F ( U ) , x ∈ R . (4)Wave trains are solutions to (4) of the form U ( x, t ) = U ∞ ( κx − ωt ) where U ∞ is 2 π -periodic in itsargument, so that κ is the spatial wave number, and ω is the temporal frequency. In the travelingcoordinate ξ = κx − ωt , wave trains are stationary solutions of U t = κ DU ξξ + ωU ξ + F ( U ) , ξ ∈ R . (5)Generically, wave trains arise as one-parameter families for which ω and κ are connected by thenonlinear dispersion relation ω = ω ∗ ( κ ) and the profile U ∞ ( ξ ; κ ) depends smoothly on κ . The groupvelocity of the wave train is defined from the nonlinear dispersion relation as c g = dωdκ : it is equal tothe speed with which perturbations are transported along the wave train in the original laboratoryframe.To prepare for our discussion of spiral waves on bounded domains, we introduce the concept ofboundary sinks which connect wave trains of (4) at x = −∞ with a Neumann boundary conditionat x = 0. We say that U ( x, t ) = U bdy ( x, ωt ) where U bdy ( x, τ ) is 2 π -periodic in τ and satisfies the5ollowing one-dimensional equation on the half line in the laboratory frame ωU τ = DU xx + F ( U ) , ( x, t ) ∈ ( −∞ , × S , (6) U x (0 , τ ) = 0 , τ ∈ S and converges to a wave train U ∞ ( κx − ωt ) with c g > x → −∞ such that (cid:12)(cid:12) U bdy ( x, · ) − U ∞ ( κx − · ) (cid:12)(cid:12) C ( S ) → . An example of a boundary sink is shown in Figure 2. (b)
Neumann Boundary t = 0t = T Space T i m e x = 0
40 1 2 3 Space x = 0 (a) Figure 2: (a) Illustration of a boundary sink. The effect of Neumann boundary conditionis highlighted by the temporal cross section on the right. (b) Comparison between planarspiral wave (left) and spiral wave on a bounded disk with Neumann boundary conditions(right). Both spirals shown on disks for comparison purposes.
We say that the reaction-diffusion system (1) has a planar spiral wave solution of the form U ( x, t ) = U ∗ ( r, φ − ωt ) where ( r, φ ) are polar coordinates if there exists an ω ∈ R and a smooth function θ ( r )with θ (cid:48) ( r ) → U ∗ satisfies (1) and | U ∗ ( r, · − ωt ) − U ∞ ( κr + θ ( r ) + · − ωt ) | C ( S ) → , as r → ∞ , where U ∞ ( κr − ωt ) is a wave train with c g >
0. Here, ω is the temporal rotational frequency and θ ( r ) acts as a phase correction to match solutions at the core with asymptotic wave trains. Thespatial wave number κ is selected by the spiral, and the wave train connects ω and κ through thenonlinear dispersion relation ω = ω ∗ ( κ ). Spiral waves are stationary solutions in the co-rotating6olar frame ( r, ψ ) = ( r, φ − ωt ) U t = D ∆ r,ψ U + ωU ψ + F ( U ) . (7)Finite domains are physically and numerically realistic, and bounded disks in particular are commoncomputational domains as they incorporate rotational symmetry properties of the spiral. Whenconsidered on B R (0), the disk of radius R centered at the origin, planar spiral waves are truncatedand solutions with positive group velocity emitted from the core are now matched with time 2 π/ω -periodic boundary sinks U bdy ( ξ, τ ). Spirals formed on B R (0) with homogeneous Neumann boundaryconditions are stationary solutions of the system in the rotating polar frame U t = D ∆ r,ψ U + ω R U ψ + F ( U ) , ( r, ψ ) ∈ [0 , R ) × S (8) U r ( R, ψ ) = 0 , ψ ∈ S . The temporal frequency of the bounded spiral converges to that on the infinite domain ω R → ω ∗ as R → ∞ . An example of spiral waves on bounded disks and the entire plane is shown inFigure 2. Note how the homogeneous Neumann boundary conditions influence the spiral near theouter boundary at the top of the spiral, similar to the boundary sink. Here, we give a review of relevant spectral and stability concepts. First, general spectral definitionsand terminology is defined. Then the spectra of the one-dimensional wave trains is described,followed by those of spirals on unbounded domains. In each case, we first consider the patterns oninfinite domains and then describe how spectra are modified by bounded domains.The spectrum Σ of a closed, densely defined linear operator L : X → X on the Banach space X isdefined as Σ = (cid:8) λ ∈ C (cid:12)(cid:12) ( L − λ ) : X → X does not have a bounded inverse (cid:9) . If L is Fredholm, the spectrum can be decomposed into the following disjoint sets [19]:Σ = Σ pt ∪ Σ Fred ∪ Σ FB where Σ pt = { λ ∈ C : L − λ is Fredholm with index 0 but not invertible } Σ FB = { λ ∈ C : L − λ is not Fredholm } Σ Fred = { λ ∈ C : L − λ is Fredholm with non-zero index } . pt is referred to as the point spectrum and contains elements called eigenvalues, which aretypically discrete. Often, the essential spectrum Σ ess is defined by Σ ess = Σ / Σ pt = Σ Fred ∪ Σ FB .The contents of each set will depend on the operator and some of these sets may be empty. Thesesets are illustrated in Figure 3. Stability of wave trains in the co-moving coordinate frame ξ = κx − ωt is analyzed by consideringthe spectrum of the operator L mv ∞ V = κ DV ξξ + ωV ξ + F U ( U ∞ ) V (9)on L ( R ) formed by linearizing (5) about the solution U ∞ . There are no non-trivial solutions V ∈ L ( R ) in the kernel of L mv ∞ − λ and the sets Σ pt and Σ Fred are empty, that is Σ pt = Σ Fred = ∅ .The spectrum Σ of L mv ∞ consists only of Fredholm borders Σ FB , which can be computed as follows.The linearization F U ( U ∞ ) is 2 π -periodic, so by Floquet theory we seek non-trivial solutions to L mv ∞ V = λV of the form V ( ξ, t ) = e νξ/κ ¯ V ( ξ ) with ¯ V ( ξ + 2 π ) = ¯ V ( ξ ) for ν ∈ C and obtain therelation L mv ∞ ( λ, ν ) ¯ V = D ( κ∂ ξ + ν ) ¯ V + ωκ ( κ∂ ξ + ν ) ¯ V + F U ( U ∞ ) ¯ V − λ ¯ V = 0 (10)which connects the temporal eigenvalues λ and spatial Floquet exponents ν . Therefore, the spec-trum Σ of L mv ∞ is given byΣ wt := (cid:8) λ ∈ C : ∃ ν ∈ i R and non-trivial 2 π -periodic ¯ V ( ξ ) so that L mv ∞ ( λ, ν ) ¯ V = 0 ∀ ξ ∈ R (cid:9) . It can be shown that Σ wt is the union of smooth curves of the form λ = λ ∞ ( ν ) with ν = iγ ∈ i R ,which are often referred to as linear dispersion curves. For each fixed λ ∈ C , equation (10) admitsfinitely many Floquet exponents ν ∈ C , and at least one Floquet exponent crosses the imaginaryaxis as λ crosses through a spectral curve.In the laboratory frame, the linearized equation is V t = DV xx + F U ( U ∞ ( κx − ωt )) V. (11)Functions of the form V ( x, t ) = e λt e νx ¯ V ( κx − ωt ) for non-trivial 2 π -periodic ¯ V ( κx − ωt ) = ¯ V ( ξ )satisfy equation (11) if and only if L lab ∞ ( λ, ν ) ¯ V = D ( κ∂ ξ + ν ) ¯ V + ω ¯ V ξ + F U ( U ∞ ( ξ )) ¯ V − λ ¯ V = 0 , (12)which defines essential spectrum curves λ = λ lab ( ν ) for ν ∈ i R . We note that the essential spectrain the co-moving (10) and laboratory frames (12) are different: comparing L mv ∞ ( λ, ν ) to L lab ∞ ( λ, ν ),8e see that the spectral curves are related via [34, 37] λ lab ( ν ) = λ ∞ ( ν ) − ωκ ν + iω(cid:96), (cid:96) ∈ Z . (13)Essential spectra computed in the laboratory frame have vertical periodic branches parameterizedby (cid:96) ∈ Z , which arise from Floquet ambiguity. Additionally, for each fixed λ ∈ C in (12), there areinfinitely many ν ∈ C and non-trivial 2 π -periodic functions ¯ V such that L lab ∞ ( λ, ν ) ¯ V = 0. We orderthe spatial eigenvalues ν j for fixed λ (cid:29) · · · ≤ Re( ν − j − ) ≤ Re( ν − j ) ≤ · · · ≤ Re( ν − ) < < Re( ν ) ≤ · · · ≤ Re( ν j ) ≤ Re( ν j +1 ) ≤ · · · (14)so that Re( ν − ) < < Re( ν ). Upon crossing an essential spectrum curve (Fredholm border), atleast one spatial eigenvalue crosses the imaginary axis. Figure 3 shows curves of essential spectrumwith insets indicating the distribution of the spatial Floquet eigenvalues ν . ⌃ abs ⌃ ess Re( )Im( ) Im( ⌫ )Re( ⌫ ) ⌃ pt Re( )Im( ) ⌃ pt ⌃ Fred ⌃ FB (b)(a) Figure 3: (a) Illustration of spectral sets for a general linear operator. In this example, theoperator
L − λ is Fredholm with index 0 in the unshaded region. Point spectrum (squares)can be found in this region and are defined for λ that make the operator not invertible.Moving to the left, the operator L − λ becomes not Fredholm on the curve Σ FB and is thenFredholm with non-zero index in the shaded region. (b) Illustration of spectra for spiralwave linear operator. Insets show distribution of spatial eigenvalues, with dots (crosses)indicating ν with initially Re( ν ) > <
0) for λ (cid:29)
1. Crossing the essential spectrum fromright to left results in one spatial eigenvalue crossing the imaginary axis and the real partbecoming positive. On absolute spectra curves, the eigenvalue that crossed the imaginaryaxis aligns with one from the initial positive set.
Stability of planar spirals can be determined similarly to the one-dimensional case by consideringthe spectrum of the operator formed by linearizing (7) about U ∗ ( r, ψ ) L ∗ V = D (cid:18) ∂ rr + 1 r ∂ r + 1 r ∂ ψψ (cid:19) V + ωV ψ + F U ( U ∗ ( r, ψ )) V (15)9nd considering the operator L ∗ on L ( R ) with domain H ( R ). Here, the spectrum Σ contains λ = 0 and λ = ± iω which arise from rotational and translational symmetries of the planar spiral.Features of Σ ess depend purely on asymptotic properties of the spiral [33, 34]. In the formal limit r → ∞ , the linear operator L ∗ becomes˜ L ∗ = D∂ rr + ω∂ ψ + F U ( U ∞ ) . (16)Eigenfunctions in the far-field limit take the form [35] V ( r, ψ ) = e νr e i(cid:96)ψ ¯ V ( κr + ψ ) , ¯ V ( ξ + 2 π ) = ¯ V ( ξ ) , (17)where radial growth or decay is characterized by the real part of the spatial eigenvalue ν and ¯ V ( ξ )is a periodic eigenfunction of the asymptotic wave train. Substitution into ˜ L ∗ V = λV gives˜ L ∗ ( λ, ν ) ¯ V = D ( κ∂ ξ + ν ) ¯ V + ω ¯ V ξ + F U ( U ∞ ) ¯ V − λ ¯ V . (18)The Fredholm borders of the essential spectrum are defined by λ = λ ∗ ( ν ) for which one spatialeigenvalue is purely imaginary [33] and we see that the far-field spiral-wave operator reduces to thecase of the laboratory frame wave train (12), that is λ ∗ ( ν ) = λ lab ( ν ). In general,Σ ess ( L ∗ ) = (cid:110) λ ∈ C : ˜ L ∗ ( λ, ν ) ¯ V = 0 has a non-trivial solution ¯ V ∈ H (cid:0) S (cid:1) for ν ∈ i R (cid:111) ∪ { λ ∈ C : L − λ ∗ is Fredholm with non-zero index } , and this set is connected to the essential spectrum of periodic wave trains in the co-moving framevia relation (13) [34]. Since ν ∈ i R , the mapping does not modify stability properties, and Σ ess ( L ∞ )and Σ ess ( L ∗ ) destabilize under the same conditions. The additional vertical periodic branches atinteger multiples of iω are distinct branches for the spiral and no longer artifacts from Floquettheory as here far-field rotational symmetry implies that if V ( r, ψ ) is an eigenfunction, then so is e i(cid:96)ψ V ( r, ψ ). As in the laboratory frame, there exists infinitely many spatial eigenvalues ν for eachtemporal eigenvalues λ , which we will order by real part, as in (14).The pertinent linear operator for spiral waves on the bounded disk B R (0) with Neumann boundaryconditions is L ∗ ,R V = D ∆ r,ψ V + ωV ψ + F U ( U ∗ ,R ) V, ( r, ψ ) ∈ [0 , R ) × S (19) V r ( R, ψ ) = 0 , ψ ∈ S . The spectrum Σ ∗ ,R of L ∗ ,R contains only point spectrum as the operator L ∗ ,R − λ is Fredholm withindex zero for all λ . Naturally, we expect that the discrete eigenvalues in Σ ∗ ,R resemble the spectraof the planar spiral wave Σ( L ∗ ), however, we will see that this is not true in general. Instead,eigenvalues in Σ ∗ ,R will converge to the union of three sets: the extended point spectrum Σ ext , the10bsolute spectrum Σ abs , and the spectrum of the boundary sink Σ bdy . We describe these sets andtheir properties below.Intuitively, the essential spectrum describes convective instabilities, in which growing perturbationsare transported away to infinity [33]. When posed on a bounded disk, convective instabilities areno longer relevant. Instead, perturbations that grow in norm at every point in space becomesignificant. These so-called absolute instabilities are captured in the limit R → ∞ by the absolutespectrum, which is defined via the far-field linear dispersion relation ˜ L ∗ ( λ, ν ) asΣ abs = { λ ∈ C : Re ν − ( λ ) = Re ν ( λ ) } . The absolute spectrum consists typically of curves that are parametrized by β = | Im ν − − Im ν | ,where β = 0 at the end points. Elements of Σ abs do not correspond to eigenvalues of L ∗ ,R butrather represent accumulation points of infinitely many discrete eigenvalues of L ∗ ,R as the domainsize R goes to infinity [33]. Note that the absolute spectrum is still defined by the limiting operator˜ L ∗ for the planar spiral wave, but it is generally distinct from the essential spectrum. The spatialeigenvalues corresponding to elements in the absolute spectrum are also illustrated in Figure 3.To explain the extended point spectrum Σ ext , we introduce spaces with exponential weight functionsin polar coordinates on R with given weight η ∈ R in the radial direction via L η (cid:0) R (cid:1) := { u ∈ L : | u | L η < ∞} , | u | L η := (cid:90) R (cid:12)(cid:12) u ( x ) e η | x | (cid:12)(cid:12) dx. For every λ / ∈ Σ abs , there exists an η ∈ R such that Re ν − ( λ ) < η < Re ν ( λ ). Using exponentiallyweighted spaces, we can then define the extended point spectrum Σ ext via [13, 33]Σ ext = (cid:8) λ ∈ C \ Σ abs : L ∗ − λ is not boundedly invertible on L η where η is such that Re ν − ( λ ) < η < Re ν ( λ ) (cid:9) . The weight permits exponential radial growth of eigenfunctions up to rate η . Exponentiallyweighted norms are equivalent on bounded domains, therefore we can expect that elements inthe extended point spectrum of planar spiral waves Σ ext ( L ∗ ) persist as eigenvalues for the operatoron each bounded disk, which we can attribute to instabilities caused by the core.Finally, the boundary conditions may contribute additional point eigenvalues, which belong tothe spectrum of the boundary sink Σ bdy defined above in equation (6). The pertinent linearizedoperator is given by L bdy V = − ωV τ + DV xx + F U ( U bdy ) V, ( x, τ ) ∈ ( −∞ , × S (20) V x (0 , τ ) = 0 , τ ∈ S . The relevant eigenvalues of the boundary sink are those that persist on finite domains [13, 33, 37].Therefore, the extended point spectrum of the boundary sink provides the eigenvalues of the spiral11aused by the boundary conditions, and we have Σ bdy = Σ ext ( L bdy ).We summarize our discussion in the following theorem from [13, Theorem 2.5.5], see [33, 38] foradditional details. Theorem 1.
The spectrum
Σ ( L ∗ ,R ) of L ∗ ,R on L ([0 , R ] × [0 , π )) with Neumann boundary con-ditions converges: Σ ( L ∗ ,R ) → Σ abs ( L ∗ ) ∪ Σ ext ( L ∗ ) ∪ Σ ext ( L bdy ) as R → ∞ . where L ∗ is the operator for a planar spiral wave on L ( R ) and L bdy is the boundary sink opera-tor. Convergence is uniform on bounded subsets of the complex plane in the symmetric Hausdorffdistance. Moreover, the multiplicity of eigenvalues in the extended point spectrum is preserved; incontrast, the number of eigenvalues, counted with multiplicity, in any fixed open neighborhood ofany point λ ∈ Σ abs ( L ∗ ) converges to infinity as R → ∞ . Therefore, on bounded domains, eigenvalues fall into one of three sets: (1) extended point spec-trum that persist under truncation, (2) eigenvalues converging to and emerging from the absolutespectrum, and (3) spectrum of the boundary sink.We note that it was proved in [36] that, under certain conditions on the asymptotic equations,isolated eigenvalues in Σ ext may emerge from absolute spectrum branch points at predictable an-gles and destabilize prior to the absolute spectrum. The location of these isolated eigenvalues ispredicted by including 1 /r curvature terms into the asymptotic problem [36, 41]. Alternans and line defects are observed on bounded domains. To investigate whether unstableeigenvalues that generate these instabilities originate from Σ abs , Σ ext , or Σ bdy , we consider spiralsformed on three domains, each of which contains eigenvalues from a portion of the spectral sets.The first domain is the standard bounded disk of radius R with homogeneous Neumann boundaryconditions, which we denote by B R (0). Here, spirals U ∗ ,R ( r, ψ ) are solutions of (8), and the spectrumof the operator L ∗ ,R V = D ∆ r,ψ V + ωV ψ + F U ( U ∗ ) V (21)on L ( B R (0)) with domain (cid:8) V ( x ) ∈ H ( B R (0)) : V r ( R, · ) = 0 (cid:9) provides information about stability.From Theorem 1, the spectrum contains contributions from the core, the far field, and the boundarysink captured by the sets Σ ext , Σ abs , and Σ bdy , respectively.Core instabilities associated with Σ ext are analyzed by computing spirals on the bounded disk radius R with non-reflecting boundary conditions. Non-reflecting boundary conditions mimic an infinite12omain by matching the spiral to the asymptotic wave train on the boundary, allowing the spiralto naturally pass through it without interference. Non-reflecting spirals U nr ( r, ψ ) are solutions to0 = D ∆ r,ψ U + ωU ψ + F ( U ) , ( r, ψ ) ∈ [0 , R ) × S U r − κU ψ , r = R, ψ ∈ S , (22)where the boundary condition is obtained by taking derivatives of the asymptotic matching condi-tion U ∗ ( r, ψ ) = U ∞ ( κr − ψ ). The linear operator for the non-reflecting spirals L R,nr is L R,nr V = D ∆ r,ψ V + ωV ψ + F U ( U nr ) V (23)which acts on eigenfunctions in (cid:8) V ( x ) ∈ H ( B R (0)) : V r ( R, ψ ) = κV ψ ( R, ψ ) (cid:9) .Finally, effects of the boundary and eigenvalues associated with Σ bdy are captured by direct compu-tation of the boundary sinks. The time 2 π/ω -periodic pattern U bdy is posed on a two-dimensionalspatiotemporal domain Ω bdy = [ − L, × S , with Neumann boundary conditions at x = 0 and2 π -periodic boundary conditions in τ . U bdy ( x, τ ) is a solution to U t τ = ω [ DU xx + F ( U )] , ( x, τ ) ∈ [ − L, × S (24) U x (0 , τ ) = 0 , τ ∈ S . (25)Stability of the boundary sink is given by considering the operator L bdy,L V = − ωV τ + DV xx + F U ( U bdy ) V (26)on the space (cid:8) V ( x, τ ) ∈ H ([ − L, × S ) : V x (0 , · ) = 0 (cid:9) . Boundary sinks for the R¨ossler andKarma models are shown in Figures 5a and 8b and will be discussed in further detail below.Boundary sinks have far-field dynamics and boundary conditions, but lack the core: conversely,non-reflecting spirals contain the core but lack outer boundary effects. On each domain, stabilityproperties are given by spectra of the operator linearized around the solution. Comparing thespectra of these three operators will indicate which region and spectral set is responsible for observedinstabilities. We expect all operators to have eigenvalues aligning along the absolute spectrum dueto the far-field dynamics, but the spectrum of L bdy,L will not contain isolated core eigenvaluesfrom Σ ext ( L ∗ ), and L R,nr will not have eigenvalues from the boundary sink. We remark that fornon-reflecting boundary conditions in L R,nr discrete eigenvalues from the far-field will still convergeto the absolute spectrum. These expectations are summarized in the following lemma [38].
Lemma 1.
The spectra of the operators defined above have the following limits, where L ∗ is thelinear operator for the planar spiral wave on L (cid:0) R (cid:1) with domain H (cid:0) R (cid:1) defined in (16), and L bdy is the boundary sink defined in (20):1. Bounded disk:
13 ( L ∗ ,R ) → Σ ext ( L ∗ ) ∪ Σ ext ( L bdy ) ∪ Σ abs ( L ∗ ) , as R → ∞
2. Non-reflecting disk:
Σ ( L R,nr ) → Σ ext ( L ∗ ) ∪ Σ abs ( L ∗ ) , as R → ∞
3. Boundary sink:
Σ ( L bdy,L ) → Σ ext ( L bdy ) ∪ Σ abs ( L ∗ ) , as L → ∞ . Numerical Methods:
The patterns and spectra of each operator are computed numerically in Matlab. Patterns areformulated as roots of equations of the form F ( U ) = 0 representing the discretized PDE posed onan appropriate domain. Solutions are found using Matlab’s built-in root finding algorithm fsolve .Periodic wave trains U ∞ ( ξ ) are found by solving0 = κ DU ξξ + ωU ξ + F ( U )on the domain ξ ∈ [0 , π ) with periodic boundary conditions U ∞ ( ξ + 2 π ) = U ∞ ( ξ ). Translationalsymmetry creates a family of solutions, and to select a unique solution and create a square systemthe phase condition 0 = (cid:90) π (cid:104) U ξ ( y ) , U old ( y ) − U ( y ) (cid:105) dy (27)and is added to F ( U ) where U old ( ξ ) is the initial guess for the wave train. One-dimensional periodicdomains are discretized using Fourier spectral differentiation matrices with N ξ = 128 grid points.Continuous spectra of the wave train are calculated through numerically continuation of the lineardispersion relation L mv ∞ ( λ, ν ) ˜ V = 0 using methods described in [29], which gives Σ FB of spiral wavesvia the relation (13).On large bounded disks, spiral waves are computed as roots of the equation0 = D (cid:18) ∂ rr + 1 r ∂ r + 1 r ∂ ψψ (cid:19) U + ωU ψ + F ( U )with appropriate boundary conditions (homogeneous Neumann or non-reflecting). The spiral an-gular frequency ω R depends on the radius R of the disk and is added as a free parameter in thespiral calculation. Rotational symmetry also creates a family of solutions, and the phase conditionfor one dimensional waves (27) is applied at r = R/ B R (0) and B nr R (0) are discretized with a fourth-order centered finitedifference scheme with N r = 200 grid points in the radial direction and periodic Fourier spectralmethods with N θ = 100 grid points in the angular coordinate. Grid sizes and discretizationsare chosen to ensure numerical accuracy and to capture a sufficient number of spiral bands for14onvergence of eigenvalues to occur, while maintaining efficient calculations. Radii of R = 125and R = 5 are used for the R¨ossler and Karma models, respectively, which results in capturing atleast three spiral bands due to the spatial wave numbers. A Neumann compatibility condition isenforced at the origin of the polar grid. As in [41], two variations of polar grids are used for thespiral and eigenvalue computations. Spiral solutions are solved on a grid of size N θ × N r , where theorigin contains N θ grid points. A grid with only one grid point at the origin is used for eigenvaluecalculations. Neumann boundary conditions on the outer radius are implemented into the finitedifference matrices via the ghost point method [40, Section 1.4].The boundary sink operator on the rectangular domain Ω bdy was discretized similarly using fourth-order centered finite differences with N s grid points in the spatial direction and a Fourier spectralmethod with N t grid points in the periodic temporal direction. Following the methods and termi-nology in [15, 22], the boundary sink is computed numerically by decomposing the domain into a“far-field” region in which the asymptotic wave train is translated in time and space and a “core”region where the Neumann boundary condition has an effect on the wave shape. Note that in thiscase the core refers to the area near the boundary. The pattern and spatial wave number of thefar-field region is fixed to match that of the spiral wave. Smooth cut-off functions of the form χ ( x ) = 1 / x − d )) match the solutions U wt ( x, τ ) in the far-field and in the core W ( x, τ ).The full solution is then given by U bdy ( x, τ ) = (1 − χ ( x )) U wt ( x, τ ) + χ ( x ) W ( x, τ ) . Substituting the form of U bdy into (6) allows us to calculate W ( x, τ ) with Newton’s method. Toaccount for numerical inaccuracies, the temporal frequency ω is set as a free parameter and anintegral phase condition is added to match the U wt and W solutions. That is, computing boundarysinks amounts to solving the system − ω∂ τ U bdy + D∂ xx U bdy + F ( U bdy ) = 0 , ( x, τ ) ∈ ( − L, × S W x (0 , τ ) = 0 , τ ∈ S (cid:90) − π/κ (cid:90) π ∂ x U wt ( x, τ ) W ( x, τ ) dτ dx = 0for W ( x, τ ) where the temporal direction is scaled to be 2 π -periodic in τ .Solutions in the far-field are obtained by translating asymptotic wave trains in time and space usingthe angular frequency and spatial wave number from the spiral and imposing the cutoff function χ ( x ). Applying (1 − χ ( x )) to the translation yields an initial condition for W ( x, τ ). Domain sizeswere selected to fit 6 periods of the wave train, which accurately captured both the Neumannboundary conditions and convergence to the far-field dynamics. Asymptotic wave trains werecomputed from the one-dimensional problem (5) using Fourier spectral methods on a periodic gridof N t points. The translation of wave train to boundary sink resulted in N s = 6 N s . To take spatialderivatives, the pattern was initially posed on a larger spatial grid of 8 periods with Neumann15oundary conditions on each end. When solving for the final pattern, the left two periods wereremoved to eliminate left-hand side boundary effects and simulate a half-infinite line. Spirals on each domain are numerically calculated for the Karma and R¨ossler models, and theinfluence of the spiral regions is determined by comparing the spectra of the three operators L ∗ ,R , L R,nr , and L bdy . At the onset of period doubling, point eigenvalues with imaginary parts approximately equal to ω + (cid:96)ω , (cid:96) ∈ Z destabilize, followed by branches of essential and then absolute spectra upon increasing µ R further. The unstable eigenfunctions are localized at the boundary (Figure 4), indicating thatline defects are a result of instabilities of the boundary conditions. The spectra of L ∗ ,R , L R,nr ,and L bdy are compared in Figure 5. As expected, all patterns have eigenvalues from the far-fielddynamics aligning along the absolute spectrum. However, only domains with boundary conditions,that is the spiral on B R (0) and boundary sink on Ω bdy contain the unstable line defect eigenvalues.Thus, the instability is confirmed to arise from the boundary conditions, and a bounded domain isnecessary for the defects to occur. 16 Re( ) I m () -0.4 -0.3 -0.2 -0.1 0 0.1 Re( ) I m () (a) (c) (b) (d) Absolute PointEssential !/ ! !/ !/ ! !/ Figure 4: R¨ossler Model: (a) Spectra for L ∗ ,R representing a stable spiral on a disk withparameter µ R = 2 and radius R = 125. Labels on right side of imaginary axis indicatehalf-multiples of angular frequency. (b) Spectra of unstable spiral, µ R = 3 .
4, (c) Spiralon bounded disk of radius R = 125 exhibiting a single stationary line defect. Parameter µ R = 3 .
4. (d) Unstable point eigenfunction responsible for line defects with µ R = 3 . λ = 0 .
043 + 0 . i = 0 .
043 + ω/ i .To further probe for influence of the boundary conditions, we can modify them by changing κ inequation (22); κ = 0 corresponds to homogeneous Neumann conditions and κ = κ ∗ is non-reflecting.Therefore, we can start at κ = 0 with the homogeneous Neumann boundary spiral from B R (0),numerically continue in κ until reaching the spatial wave number of the spiral κ ∗ and track theevolution of an unstable eigenvalue. The spiral is already formulated as a root finding problem, andthe eigenvalue problem can be as well by solving L R,nr V − λV = 0 . Each continuation stage is a2-step process. First, a new spiral with updated boundary conditions is computed, and second thelinearized operator L R,nr is modified and an eigenpair ( λ, V ) is computed. Starting the continuationat κ = 0 allows the unstable eigenfunction from L ∗ ,R to be used in the first continuation step. Ifthe eigenvalue is unchanged with the boundary continuation, then it does not originate from theboundary sink. 17 a) N e u m a nn B o und a r y t = 0t = T ⇠ = 0 Space T i m e Neumann Spiral ( ) + Non-Reflecting Spiral( ) Boundary Sink ( ) L ⇤ ,R L ⇤ , nr L bdry (b) -0.05 0 0.05 Re( ) -2-1012 I m () Figure 5: R¨ossler Model: (a) Image of boundary sink. Neumann boundary conditions onthe right at ξ = 0. Domain is periodic in time (vertical) direction. (b) Spectra of operators L ∗ ,R , L R, nr , L bdy .Figure 6 shows the evolution of the point eigenvalue during the κ -continuation. The eigenvaluechanges with κ , first tracking along the essential spectrum, and then jumping on the absolutespectrum. Increasing κ corresponds to a mixed boundary condition and results in different shapesof unstable eigenfunctions, demonstrating that the boundary conditions will change the observedinstability. The eigenfunctions in Figure 5(b) show the transition from localization at the boundaryto localization at the core as κ is increased from 0 to κ ∗ . Similar results are obtained for unstablepoint eigenvalues at other multiples of i ω + iω(cid:96) as these eigenvalues arise due to the asymptotic e i(cid:96)ψ V ( r, ψ ) symmetry of the eigenfunctions. 18 ! * ! -0.06 -0.04 -0.02 0 0.02 0.04 0.06 Re( ) I m () !/ ! !/ (a)(b) NeumannNon-Reflecting
Figure 6: R¨ossler Model: (a) Boundary condition continuation of line defect eigenvalue(black curve), starting with κ = 0 on the right and ending on absolute spectrum with κ = κ ∗ . (b) Eigenfunctions from continuation. Locations of eigenvalues indicated by greencircles on eigenvalue continuation in (a). Eigenfunctions are on a disk of radius R = 125. As the bifurcation parameter µ K is increased above one in the Karma model, the essential spectrumdestabilizes in an Eckhaus instability, followed by a single complex-conjugate pair of eigenvalueswith imaginary part near 3 ω/
2. Meandering and alternans appear with the Hopf bifurcation fromthe point eigenvalues. Shown in Figure 7d, the unstable point eigenfunction, and hence the form ofthe instability, has highest magnitude at the boundary of the spiral bands, leading to the observedalternans.The single pair of eigenvalues suggests they arise from Σ ext and are instabilities of the core. Inthis case, comparison of the three operators indicates alternans eigenvalues are present in thespiral spectra for L ∗ ,R and L R,nr , but are absent in the boundary sink L bdy . Modification of theboundary conditions in B R,nr (0) by continuing κ results in no change to the alternans eigenvalueor eigenfunction. Furthermore, the pair is not emitted from the absolute spectrum; continuationof the absolute spectrum branch point λ bp and alternans eigenvalue λ A in parameter µ K shows thedifference Re( λ bp ) − Re( λ A ) is positive over an appropriate range of parameter values (Figure 8c).Thus, the point eigenvalues causing the alternans instability stem instead from the unstable pairof eigenvalues originating from Σ ext affiliated with the core.19 a) (d) !/ ! !/ (c) Time !/ ! !/ (b) Absolute PointEssential
Figure 7: Karma Model: (a) Spectra for stable spiral, µ K = 0 .
6. Labels on right side ofimaginary axis indicate half-multiples of angular frequency. (b) Spectra of unstable spiral, µ K = 1 .
4, (c) Development of alternans in time evolution of spiral on bounded square withhomogeneous-Neumann boundary conditions. Parameter µ K = 1 .
4. Square of side length16cm. (d) Unstable point eigenfunction responsible for alternans. Corresponds to eigenvalue λ = 2 . . i ≈ . ω/ i , µ K = 1 .
4. Domain radius R = 5 with homogeneous-Neumannboundary conditions. More can be said about the alternans eigenfunction. To leading order, spiral eigenfunctions are ofthe form (17), and as the unstable alternans point eigenvalue passes through the essential spectrum,the eigenfunction inherits properties of the continuous spectrum [35]. In particular, Re ν is smallwhen the eigenvalue is near the essential spectrum, with Re ν positive (negative) for the eigenvalueto the left (right) of the continuous spectrum curve. Small Re ν indicates little radial growth, andIm ν = γ is set by the essential spectrum point. The wave train eigenfunction is determine by theeigenfunction on the essential spectrum, V ess . Using these observations, the spiral eigenfunction toleading order in the radius is given by V ( r, ψ ) = e iγr e (cid:96)ψ V ess ( κr + ψ ) . (28)A numerically constructed eigenfunction is shown at the top of Figure 9a. Note that the derivativeof the wave train ∂ ξ U ∞ ( ξ ) is the eigenfunction on the essential spectrum with eigenvalue λ = i(cid:96)ω , (cid:96) ∈ Z . The alternans eigenfunction crosses Σ ess near one of these points, and therefore V ess ( ξ ) isclose to U (cid:48)∞ ( ξ ) (Figure 9b). The derivative of the wave train is the highest at the wave frontsand backs and it is this shape that leads to the changing width of the spiral bands and form of20 .2 0.4 0.6 0.8 1 1.2 R R ea l ( A - bp ) (a) (c) t = 0 N e u m a nn B o und a r y t = T ⇠ = 0 T i m e Space (b)
Neumann Spiral ( ) + Non-Reflecting Spiral( ) Boundary Sink ( ) L ⇤ ,R L ⇤ , nr L bdry -6 -4 -2 0 2 4 Re( ) -100-50050100 I m () Figure 8: Karma Model: (a) Spectra of operators L ∗ ,R , L R, nr , L bdy . (b) Image of boundarysink. Neumann boundary conditions on right at ξ = 0. Periodic in time (vertical) direction.(c) Distance between real part of alternans point eigenvalue and branch point of absolutespectrum as a function of parameter µ R .alternans. Moreover, the structure of the constructed eigenfunction is in good agreement with thealternans eigenfunction from L ∗ ,R which is reproduced on the bottom of Figure 9a.When the alternans eigenvalue is to the left of the essential spectrum, approximately µ K < .
4, theoverall shape of the eigenfunction is comparable to those shown in Figure 9a, but there is slightradial growth toward the boundary, corresponding to a spatial eigenvalue, ν , with a small positivereal part. Radial growth of the alternans eigenfunction for several values of bifurcation parameter µ K is visible in Figure 9c. As the eigenvalue approaches the essential spectrum, radial growthdecreases and is near zero for µ K near 1.4, as seen in the solid black curve in Figure 9c. Alternansappear when this eigenvalue first destabilizes, but the strength of the instability on the core andfar-field regions is determined by the proximity of the point eigenvalue to the essential spectrum.21 Radius m a x ( | V | ) u ⌃ ess eigenfunction (a) (c)(b) Figure 9: (a) Karma model alternans eigenfunction constructed from essential spectrumdata. (b) Comparison between derivate of wave train (blue) and essential spectrum eigen-function (orange dots), (c) Radial growth of alternans eigenfunctions. Legend indicatesvalue of bifurcation parameter µ K . Instabilities in patterns observed on finite domains may be initiated by unstable eigenvalues stem-ming from a variety of sources. Determining where unstable eigenvalues originate from yields insightinto what creates the accompanying instabilities and what their spatial shape looks like. Here, wepresent a methodology for unfolding the origin of these eigenvalues and apply it to the specific caseof period-doubling instabilities of spiral waves on bounded disks. The technique of comparing thespectra of the three operators can be applied more generally to pattern forming systems on anydomain, as long as the patterns of interest can be computed as roots of an appropriate system.Our results predict that line defects in the R¨ossler model will only be seen on bounded domains andthat the shape and type of boundary conditions will likely affect the structure of the instability.Furthermore, the interaction of the outer bands of multiple spirals can induce a non-trivial boundarycondition between the spirals, and line defects or similar structures may be generated in thesesituations. Therefore, for instabilities of the boundary sink, an accurate representation of theboundary conditions and domain is a necessary factor when matching models and theory withexperiments.We find that alternans are a product of unstable eigenvalues in the extended point spectrum22ssociated with the spiral core, implying that, as long as the core does not directly interact withthe boundary, the shape of the domain and the precise form of boundary conditions are insignificantfactors in the spiral stability and formation of alternans. This result has direct impacts for cardiacdynamics in that conclusions for the development of alternans on bounded disks can be extended toirregular and complex geometries such as the heart. Unstable alternans eigenfunctions do exhibitslight radial growth or decay depending on the value of parameter µ K , which may influence whichregions of the spiral are impacted the most by the instability or how prevalent alternans are on asmall domain. Our results are consistent with [23, 24] which finds that alternans development ismost sensitive to perturbations near the core. Both results also suggest that a perturbation mustbe applied to the core region to destabilize an alternans spiral.Despite the bounded domain, the essential spectrum provides useful information about the sourceof the alternans instability. Alternans on a ring were previously attributed to destabilization ofthe continuous spectra [3]. We find that, on a bounded disk, an unstable point eigenvalue pass-ing through the Eckhaus unstable essential spectrum is fundamental to the alternans structure.Translational symmetry of the wave train ensures existence of the eigenvalue-eigenfunction pair( λ = 0 , V = U (cid:48)∞ ), meaning that the essential spectrum curve that passes through the origin for oneparameter set must contain the origin for all parameter values. Therefore, these branches generi-cally destabilize through an Eckhaus instability. Point eigenvalues that interact with these essentialspectrum branches acquire an eigenfunction with shape close to U (cid:48)∞ which impacts the wave frontsand backs leading to the observed form of the alternans instability. The unstable essential spectrumalso implies that the associated wave trains on R undergo an instability as well.In the R¨ossler model, point eigenvalues do not interact with unstable essential spectra branchesand unstable continuous spectra branches do not pass through the origin. Our results suggestthat alternans instability will appear if a point eigenvalue crosses essential spectrum that hasdestabilized through an Eckhaus instability. Future work includes investigation into whether anEckhaus instability is necessary or sufficient for the formation of alternans. If the continuousspectrum can be attributed to formation of alternans, the 1D computation of the essential spectrumprovides a tractable tool for analysis of more realistic and complex models.There are a number of differences between the spectra of the two models. In the R¨ossler system,countably many discrete eigenvalues destabilize ahead of distinct period-doubling essential and ab-solute spectrum branches. In contrast, it is a single pair of unstable complex conjugate eigenvalueswith imaginary part approximately 3 ω/ i ω + iω(cid:96) , (cid:96) ∈ Z which may take the form of smooth curvesthat intersect with or coincide with symmetry lines, or disjoint branches that do not intersect [37].The first case is observed in the R¨ossler system, and latter in Karma. Furthermore, in the R¨osslersystem, the point Im( λ ) = ω + ω(cid:96) on the continuous spectra has spatial eigenvalue Im( ν ) = κ/ ω/ (cid:96)ω . On the other23and, the disjoint absolute spectrum branches in the Karma model result in an unstable absolutespectrum contributing many frequencies close to, but not specifically period-doubling.In real cardiac systems, alternans lead to spiral break up, providing evidence they originate througha subcritical bifurcation [14, 31]. Numerical studies are inconclusive and show both immediatebreak up and short time alternans persistence [9, 20, 21]. In [17], alternans are analytically shownto originate in a subcritical Hopf bifurcation, but this analysis is limited as it relies on a specificnormal form for systems near a saddle node of a traveling wave which is not satisfied in all alternansgenerating models. The debate of a sub- versus super-critical bifurcation can be investigated inmodels by determining whether an alternans spiral is stable when considered as a time-periodicthree-dimensional structure on a bounded disk. Acknowledgements.
Dodson was supported by the NSF through grants DMS-1148284 and1644760. Sandstede was partially supported by the NSF through grant DMS-1714429.
The standard form of the Karma model in the laboratory frame for x ∈ R is E t = γ ∆ E + 1 τ E (cid:18) − E + (cid:0) E ∗ − n M (cid:1) (1 − tanh( E − E h )) E (cid:19) n t = δ ∆ n + 1 τ n (cid:18) − e − Re θ s ( E − E n ) − n (cid:19) , where E = E ( x, t ) represents membrane voltage and n = n ( x, t ) takes the place of a slower gatingvariable. In the notation of the main paper, the variables u and v represent E and n , respectively.The Heaviside function has been replaced by the smoothed function θ s ( u ) = (1 + tanh( su )) /
2. Fullparameter values are given in Table 1. The bifurcation parameter (typically called the restitutionparameter), µ K = Re is increased from 0.6 to 1.4 and controls recovery properties of the excitablemedia. All other parameters are held fixed in our study.In the laboratory frame, the R¨ossler model is given by u t = δ ∆ u − v − wv t = δ ∆ v + u + avw t = δ ∆ w + uw − cw + b. The bifurcation parameter c is increased from 2 to 3.4, with line defects appearing as µ R = c passesthrough 3. Parameters are given in Table 1. 24arma R¨ossler γ = 1 . δ = 0 . δ = 0 . δ = 0 . τ E = 0 . δ = 0 . τ n = 0 . a = 0 . E ∗ = 1 . b = 0 . M = 4 µ R = c ∈ [2 , . s = 4 ω ∈ [1 . , . E h = 3 E n = 1 µ K = Re ∈ [0 . , . ω ∈ [60 . , . ω is selected by spirals and given intervalsalign with bifurcation parameter. References [1] D. Allexandre and N. F. Otani. Preventing alternans-induced spiral wave breakup in cardiactissue: An ion-channel-based approach.
Physical Review E , 70(6):061903–16, Dec. 2004.[2] S. Alonso, M. B¨ar, and B. Echebarria. Nonlinear physics of electrical wave propagation in theheart: a review.
Reports on Progress in Physics , 79(9):096601–57, Aug. 2016.[3] M. B¨ar and L. Brusch. Breakup of spiral waves caused by radial dynamics: Eckhaus and finitewavenumber instabilities.
New Journal of Physics , 6:5–5, Jan. 2004.[4] M. B¨ar and M. Or-Guil. Alternative Scenarios of Spiral Breakup in a Reaction-DiffusionModel with Excitable and Oscillatory Dynamics.
Physical Review Letters , 82(6):1160–1163,Feb. 1999.[5] D. Barkley. Euclidean symmetry and the dynamics of rotating spiral waves.
Physical ReviewLetters , 72(1):164–167, Jan. 1994.[6] S. Bauer, G. R¨oder, and M. B¨ar. Alternans and the influence of ionic channel modifications:Cardiac three–dimensional simulations and one-dimensional numerical bifurcation analysis.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 17(1):015104–16, Mar. 2007.[7] E. M. Cherry and F. H. Fenton. Suppression of alternans and conduction blocks despite steepAPD restitution: electrotonic, memory, and conduction velocity restitution effects.
AJP: Heartand Circulatory Physiology , 286(6):H2332–H2341, June 2004.[8] E. M. Cherry, F. H. Fenton, T. Krogh-Madsen, S. Luther, and U. Parlitz. Introduction toFocus Issue: Complex Cardiac Dynamics.
Chaos: An Interdisciplinary Journal of NonlinearScience , 27(9):093701–9, Sept. 2017. 259] M. Courtemanche, L. Glass, and J. P. Keener. Instabilities of a propagating pulse in a ring ofexcitable media.
Physical Review Letters , 70(14):2182–2185, Apr. 1993.[10] E. Cytrynbaum and J. P. Keener. Stability conditions for the traveling pulse: Modifying therestitution hypothesis.
Chaos: An Interdisciplinary Journal of Nonlinear Science , 12(3):788–799, Sept. 2002.[11] B. Echebarria and A. Karma. Spatiotemporal control of cardiac alternans.
Chaos: An Inter-disciplinary Journal of Nonlinear Science , 12(3):923–930, Sept. 2002.[12] B. Echebarria and A. Karma. Amplitude equation approach to spatiotemporal dynamics ofcardiac alternans.
Physical Review E , 76(5):051911–23, Nov. 2007.[13] B. Fiedler and A. Scheel. Spatio-Temporal Dynamics of Reaction-Diffusion Patterns. In
Trendsin Nonlinear Analysis , pages 23–152. Springer, 2000.[14] L. H. Frame and M. B. Simson. Oscillations of conduction, action potential duration, and re-fractoriness. A mechanism for spontaneous termination of reentrant tachycardias.
Circulation ,78(5):1277–1287, Nov. 1988.[15] R. Goh and A. Scheel. Pattern-forming fronts in a Swift-Hohenberg equation with directionquenching - parallel and oblique stripes.
Journal of the London Mathematical Society , 98:104–128, Mar. 2018.[16] A. Goryachev, H. Chat´e, and R. Kapral. Synchronization Defects and Broken Symmetry inSpiral Waves.
Physical Review Letters , 80(4):873–876, Jan. 1998.[17] G. A. Gottwald. Bifurcation analysis of a normal form for excitable media: Are stable dynam-ical alternans on a ring possible?
Chaos: An Interdisciplinary Journal of Nonlinear Science ,18(1):013129–17, Mar. 2008.[18] M. R. Guevara, G. Ward, A. Shrier, and L. Glass. Electrical alternans and period-doublingbifurcations.
IEEE Computers in Cardiology , pages 1–5, Sept. 1984.[19] T. Kapitula and K. Promislow. Essential and Absolute Spectra. In
Spectral and DynamicalStability of Nonlinear Waves , pages 39–74. Springer, 2013.[20] A. Karma. Spiral Breakup in Model Equations of Action Potential Propagation in CardiacTissue.
Physical Review Letters , 71(7):1103–1106, Aug. 1993.[21] A. Karma. Electrical alternans and spiral wave breakup in cardiac tissue.
Chaos: An Inter-disciplinary Journal of Nonlinear Science , 4(3):461–472, 1994.[22] D. J. B. Lloyd and A. Scheel. Continuation and Bifurcation of Grain Boundaries in the Swift–Hohenberg Equation.
SIAM Journal on Applied Dynamical Systems , 16(1):252–293, Jan. 2017.2623] C. D. Marcotte and R. O. Grigoriev. Unstable spiral waves and local Euclidean symme-try in a model of cardiac tissue.
Chaos: An Interdisciplinary Journal of Nonlinear Science ,25(6):063116–17, June 2015.[24] C. D. Marcotte and R. O. Grigoriev. Adjoint eigenfunctions of temporally recurrent single-spiral solutions in a simple model of atrial fibrillation.
Chaos: An Interdisciplinary Journal ofNonlinear Science , 26(9):093107–15, Sept. 2016.[25] D. Noble. A Modification of the Hodgkin-Huxley Equations Applicable to Purkinje FibreAction and Pace-Maker Potentials.
Journal of Physiology , 160:317–352, 1962.[26] J.-S. Park and K. J. Lee. Line-defects-mediated complex-oscillatory spiral waves in a chemicalsystem.
Physical Review E , 73(6):535–15, June 2006.[27] J. M. Pastore, S. D. Girouard, K. R. Laurita, F. G. Akar, and D. S. Rosenbaum. MechanismLinking T-Wave Alternans to the Genesis of Cardiac Fibrillation.
Circulation , 99(10):1385–1394, Mar. 1999.[28] V. Perez-Munuzuri, R. Aliev, B. Vasiev, V. Perez-Villar, and V. I. Krinsky. Super-spiralstructures in an excitable medium.
Nature , 353(6346):740–742, Oct. 1991.[29] J. D. M. Rademacher, B. Sandstede, and A. Scheel. Computing absolute and essential spectrausing continuation.
Physica D: Nonlinear Phenomena , 229(2):166–183, 2007.[30] A. D. J. Robertson and J. F. Grutsch. Aggregation in Dictyostelium discoideum.
Cell , 24:603–611, June 1981.[31] D. S. Rosenbaum, L. E. Jackson, J. M. Smith, H. Garan, J. N. Ruskin, and R. J. Cohen.Electrical Alternans and Vulnerability to Ventricular Arrhythmias.
New England Journal ofMedicine , 330(4):235–241, Jan. 1994.[32] B. Sandstede. Stability of Travelling Waves. In B. Fiedler, editor,
Handbook of DynamicalSystems, Volume 2 , pages 983–1055. Elsevier, 2002.[33] B. Sandstede and A. Scheel. Absolute and convective instabilities of waves on unbounded andlarge bounded domains.
Physica D: Nonlinear Phenomena , pages 233–277, 2000.[34] B. Sandstede and A. Scheel. Absolute versus convective instability of spiral waves.
PhysicalReview E , pages 1–7, 2000.[35] B. Sandstede and A. Scheel. Superspiral Structures of Meandering and Drifting Spiral Waves.
Physical Review Letters , 86(1):171–174, 2001.[36] B. Sandstede and A. Scheel. Curvature effects on spiral spectra: Generation of point eigenval-ues near branch points.
Physical Review E , 73(1):016217–8, Jan. 2006.2737] B. Sandstede and A. Scheel. Period-Doubling of Spiral Waves and Defects.
SIAM Journal onApplied Dynamical Systems , 6(2):494–547, 2007.[38] B. Sandstede and A. Scheel. Spiral waves: linear and nonlinear theory.
In preparation , 2019.[39] B. Sandstede, A. Scheel, and C. Wulff. Bifurcations and Dynamics of Spiral Waves.
Journalof Nonlinear Science , pages 1–40, June 1999.[40] J. W. Thomas.
Numerical Partial Differential Equations: Finite Difference Methods , volume 22of
Texts in Applied Mathematics . Springer New York, New York, NY, 1995.[41] P. Wheeler and D. Barkley. Computation of Spiral Spectra.
SIAM Journal on Applied Dy-namical Systems , 5(1):157–177, Jan. 2006.[42] N. Wiener and A. Rosenblueth. The mathematical formulation of the problem of induction ofimpulses in a network of connected excitable elements, specifically in cardiac muscle.
Archivsdel instituto de Cardiologia de Mexico , 16(3):205–265, 1946.[43] A. Winfree. Electrical turbulence in three-dimensional heart muscle.
Science , 266(5187):1003–1006, Nov. 1994.[44] A. T. Winfree. Spiral Waves of Chemical Activity.
Science , 175:634–636, 1972.[45] M. Yoneyama, A. Fujii, and S. Maeda. Wavelength-Doubled Spiral Fragments in PhotosensitiveMonolayers.
Journal of the American Chemical Society , 117(31):8188–8191, Aug. 1995.[46] A. N. Zaikin and A. M. Zhabotinsky. Concentration Wave Propagation in Two-dimensionalLiquid-phase Self-oscillating System.