Memory effects, transient growth, and wave breakup in a model of paced atrium
MMemory effects, transient growth, and wave breakup in a model of pacedatrium
Alejandro Garz´on
1, 2, a) and Roman O. Grigoriev Department of Mathematics, Universidad Sergio Arboleda, Bogot´a 110221, Colombia School of Physics, Georgia Institute of Technology, Atlanta, Georgia 30332-0430,USA (Dated: 25 September 2018)
The mechanisms underlying cardiac fibrillation have been investigated for over a century, but we are stillfinding surprising results that change our view of this phenomenon. The present study focuses on the transitionfrom normal rhythm to atrial fibrillation associated with a gradual increase in the pacing rate. While someof our findings are consistent with existing experimental, numerical, and theoretical studies of this problem,one result appears to contradict the accepted picture. Specifically we show that, in a two-dimensional modelof paced homogeneous atrial tissue, transition from discordant alternans to conduction block, wave breakup,reentry, and spiral wave chaos is associated with transient growth of finite amplitude disturbances rather thana conventional instability. It is mathematically very similar to subcritical, or bypass, transition from laminarfluid flow to turbulence, which allows many of the tools developed in the context of fluid turbulence to beused for improving our understanding of cardiac arrhythmias.
Atrial fibrillation is a common type of cardiacarrhythmia characterized by spatially uncoordi-nated high-frequency patterns of electrical exci-tation waves. Several different explanations ofhow fibrillation arises from a highly coordinatedand regular normal rhythm have been offered.The most common is a dynamical scenario thatinvolves several stages. First, a milder type ofarrhythmia known as alternans, characterized byregular spatial and temporal modulation in theduration of the excitation, is created as a resultof a period doubling or Hopf bifurcation. In thesecond stage, the modulation grows until someportion of an excitation wave fails to propagate(which is known as conduction block), producingwave breakup. In later stages, spiral waves seg-ments, or wavelets, form and undergo subsequentbreakups until a complicated dynamical equilib-rium is established, with wavelets continually an-nihilating and reemerging. The current paradigmis that each of these stages represents a separatebifurcation that leads to either a change in thestability, or the disappearance, of the solutiondescribing a previous stage. Our analysis of asimple model of paced cardiac tissue shows thattransitions between different stages can happenwithout any bifurcations taking place. In partic-ular, conduction block can be caused by strongtransient amplification of finite amplitude distur-bances, which is a flip side of a very well-knownmemory effect that describes the sensitivity of theasymptotic dynamics of paced cardiac tissue tothe entire pacing history. a) Electronic mail: [email protected]
I. INTRODUCTION
Fibrillation, a state of uncoordinated contraction ofthe heart muscle, can occur in the ventricles, leading tosudden cardiac arrest, a major cause of death , or theatria, where it creates significant health complications .Understanding of the transition from the normal car-diac rhythm, also known as the 1:1 response, to fib-rillation is critical for the development of therapies forthe prevention or termination of this dangerous arrhyth-mia. Abundant experimental evidence points to the roleof alternans , a beat-to-beat alternation in the actionpotential duration (APD) at high excitation rates, alsoknown as the 2:2 response, as a precursor of fibrillationboth in the ventricles and in the atria . Strong cor-relation between alternans and fibrillation was also foundin the simulations employing high-fidelity models of two-dimensional ventricular and atrial tissue.In tissue, alternans manifests as spatial modulation inthe width of the excitation waves , which leads to dy-namical variation in the tissue refractoriness. In partic-ular, spatially discordant alternans, in which APD al-ternates out of phase in different regions of the heart,markedly enhances dispersion of refractoriness, increas-ing the likelihood of conduction block and, in two- orthree-dimensional tissue, reentry . Theoretical analysisof propagating excitation waves in one dimension basedon amplitude equations show that voltage-driven al-ternans are described by stable periodic or quasiperiodicsolutions, with the tissue size and dispersion of conduc-tion velocity determining whether the alternans is con-cordant or discordant. The same is true of calcium-drivenalternans , which qualitatively behave like their voltage-driven cousins. Theoretical studies of alternans in higherdimensions are too limited in scope to make any spe-cific predictions regarding the transition from alternansto fibrillation, but based on experiments and numericalsimulations we know that conduction block that leads a r X i v : . [ q - b i o . T O ] A p r to wave breakup plays a key role.As the pacing rate or the tissue size in increased, con-duction block is predicted to arise when alternans be-comes unstable leading to unbounded growth of the am-plitude of modulation . Conduction block leads toa time-periodic response such as 2:1 (every other exci-tation wave is blocked) in one dimension, but typicallygenerates wave breakup and transition to spiral wavechaos (SWC) in two dimensions or scroll wave chaosin three dimensions (respectively, atrial and ventricularfibrillation) . As we show in the present study, there isanother possibility: conduction block can be triggered byvery small, but finite, perturbations of stable discordantalternans.Although this result appears counterintuitive at first,it is not completely unexpected. It is well-known thatpaced cardiac tissue displays multistability in zero ,one , and two dimensions. Transitions between dif-ferent stable attractors in spatially extended dissipativenonequilibrium systems are in fact not uncommon. Tran-sition from stable laminar flow to turbulence in shearfluid flows is perhaps the oldest and best known exam-ple. In the case of fluid turbulence, this is known as abypass transition and involves strong transient ampli-fication of small disturbances .Although no direct evidence of the role of transient am-plification in the dynamics of cardiac tissue is availableat present, indirect evidence is provided by the so-calledmemory effect . This term refers to the dependenceof asymptotic states of paced tissue preparations on theentire pacing history rather than just the final pacing in-terval. For instance, in canine ventricles, a pacing-downprotocol (with pacing interval decreased in steps of 50ms to 10 ms) produced a different stationary pattern ofalternans than a protocol with constant pacing rate ap-plied to quiescent tissue, for the same asymptotic rate .Similar dependence was found in normal rabbit hearts :while decreasing the pacing interval by 10 ms steps failedto induce discordant alternans, this state was reached byusing smaller steps of 5 ms to 2 ms.The objective of this paper is to establish a connec-tion between transient amplification, memory effects,and the transition to atrial fibrillation in a simple two-dimensional model of cardiac tissue featuring discordantalternans. This paper is structured as follows. SectionII describes the model used in this study. The resultsand their discussion are presented in Sections III and IV,respectively. II. MODEL OF PACED ATRIAL TISSUE
Since our focus here is on fundamental mechanismsrather than detailed comparison with experiments, wechose a simple monodomain model of atrial tissue. Ithas the form of a reaction-diffusion partial differentialequation (PDE) ∂ t u = D ∇ u + f ( u ) + I p , (1) where the field u = [ u, v ]( r , t ) describes the state of car-diac cells (cardiomyocytes), D is a diagonal matrix of dif-fusion coefficients that represents the coupling betweenadjacent cells, f ( u ) is the ionic model that describes thedynamics of individual cells, and I p is the (scaled) densityof the external pacing current.Furthermore, we chose the ionic model introduced byKarma , which captures the essential alternans instabil-ity responsible for initiating and maintaining complex ar-rhythmic behaviors. To make it differentiable, the modelhas been modified such that f ( u ) = (cid:20) ( u ∗ − v M ) { − tanh( u − } u / − u(cid:15) { β Θ α ( u −
1) + Θ α ( v − v − − v } (cid:21) , (2)where Θ α ( u ) = [1 + tanh( αu )] / u is a scaled trans-membrane voltage, and v is a gating variable. Un-less otherwise stated, the parameters values used were u ∗ = 1 . M = 4, (cid:15) = 0 . D = 1 . × − cm / ms, D = D / α = 32, and β = [1 − exp( − R )] − , where R = 1 .
273 is the restitution parameter. We also used no-flux boundary conditions ∇ u · ˆ n = 0, where ˆ n is a unitvector normal to the domain boundary, as a physiologi-cally accurate representation of the boundary of excitabletissue.At the chosen value of the restitution parameter, in twodimensions, this model reproduces the transition fromnormal rhythm to alternans to fibrillation observed inexperiments as pacing frequency is increased. Yet, withjust two variables, as opposed to several tens of variablesin the most complex models , the Karma model is sim-ple enough to facilitate the computationally demandinganalysis presented in this work.All numerical simulations were carried out on a squaredomain of side length L = 2 .
489 cm, which is close tothe minimal size required for sustained SWC in Karmamodel with the present choice of parameters. The PDE(1) was discretized in space using finite differences on a96 ×
96 computational grid, which corresponds to the sizeof one computational cell ∆ x = ∆ y = 0 . u and v on this grid define the state of the sys-tem as a point in a 18432-dimensional state space. Thediscretized equations were solved using 4th order Runge-Kutta method with time step ∆ t = 0 .
01 ms chosen to bemuch smaller than the shortest characteristic time scaleof the model, τ u = 2 . × p ( t ) = (cid:26) I , ≤ t ≤ , , otherwise . (3)Denoting t n the time at which the n -th pulse is applied, n T n ( m s ) FIG. 1. Pacing protocol A (solid line) and values at whichpacing protocol B differs from A (dashed line). the pacing current can be written as I p ( r , t ) = ∞ (cid:88) n =0 p ( t − t n ) g ( r ) (cid:20) (cid:21) , (4)where g ( r ) = 1 for paced grid points and zero otherwise. III. RESULTSA. Wave breakup induced by fast pacing
We investigated stability of the excitation waves pro-duced by pacing with gradually decreasing pacing inter-val T n = t n − t n − , n = 1 , , ... , starting at t = t = 0.For some tissue models, a decrease in the pacing inter-val does not produce wave breakup unless (i) the tissueis heterogeneous or (ii) an external current is applied ata location other than the pacing site (ectopic beat) .However, for the Karma model considered here, neitherof these conditions is necessary for wave breakup.Let us define T [ i,j ] , j > i , to be a sequence of pacingstimuli with a constant interval: T i = T i +1 = ... = T j .Consider the pacing protocol A: T [1 , = 120 ms, fol-lowed by T [11 , = 110 ms, T [21 , = 100 ms, T [31 , =90 ms, and T [41 , = 80 ms (cf. Fig. 1) with initial con-dition corresponding to the stable uniform equilibrium u = [0 , T [41 , = 85 ms (dashed line in Fig. 1). Forprotocol B, each pacing stimulus generates a travelingwave with nearly circular wavefront and waveback (cf.Fig. 2(a)) that is absorbed at the bottom and right do-main boundaries. The asymptotic state is time-periodic,with period 2 T = 160 ms, and corresponds to discordantalternans as we will see below. During the first 55 pacingintervals, pacing protocol A produces qualitatively thesame dynamics as protocol B. On the 56th pacing inter-val, however, conduction block causes wave breakup (cf.Fig. 2(b)), creating two spiral waves on the 57th inter-val (cf. Fig. 2(c)), which undergo additional breakupsand eventually transition to SWC (cf. Fig. 2(d)). SWCcorresponds to atrial fibrillation; both feature re-entrantwaves and are sustained irrespective of the pacing.So how can we understand the difference in the out-comes of the two protocols? Although this difference is a (a) (b)(c) (d) FIG. 2. Voltage component, u ( r , t ), of excitation wavesduring protocol A at (a) t = t + 68 ms, (b) t = t + 72ms, (c) t = t + 28 ms, and (d) t = t + 32 ms. The initialcondition at t = 0 is the rest state. The pacing site is locatedat the upper-right corner. The same color bar is used in allsubsequent figures showing the voltage field. clear indication that the memory effect plays an impor-tant role in our system, this does not explain much. Thememory is actually quite long: indeed, 15(!) pacing stim-uli preceding conduction block in protocol A were deliv-ered at the same asymptotic rate. Given that the asymp-totic alternans state is stable , it would have been naturalto expect that any deviation present at t = t would allbut disappear after such a long period. In fact, as weshow below, infinitesimal initial perturbations about thisstable alternans state can grow , albeit transiently, by sev-eral orders of magnitude over long time intervals. Small,but finite, initial perturbations can grow enough to causea transition to a completely different (here chaotic) at-tractor, as it happens for protocol A . This is a signatureof multistability often associated with subcritical insta-bilities; similar phenomenology is found during bypasstransition to turbulence in fluid flows . B. Periodic solutions
Let us start the analysis of the problem with a descrip-tion of the various solutions to (1). To explicitly expressthe dependence of a solution u ( r , t ) of (1) on the initialcondition u ( r ), we define the time evolution operator u ( r , t ) = Φ t [ u ( r )] . (5) FIG. 3. Voltage component of the periodic solutions repre-senting the 1:1 response (left column), 2:2 response (middlecolumn), and 2:1 response (right column). These are shown attimes t = ( n + 0 . T, ( n + 0 . T, ( n + 1 . T, and ( n + 1 . T ,where n = 1 , , , ... (from top to bottom row, respectively). For constant pacing, T [1 , ∞ ] = T , various time-periodicsolutions can be defined as fixed points or k -cycles ofthe Poincar´e map u n ≡ u ( r , nT ) = Φ T [ u n − ], which de-scribes the discrete-time dynamics on the Poincar´e sec-tion t = nT . States u k ( r , t ) with different periods corre-sponds to solutions of the nonlinear equationΦ kT [ u ] = u (6)with different k . In particular, the 1:1 response u ( r , t )corresponds to the fixed point of the Poincar´e map u = Φ T [ u ] = u ( r , k = 2 appear. The 2-cycle { u , u } , where u = Φ T [ u ] = u ( r ,
0) and u = Φ T [ u ] = u ( r , T ), corresponds to the alter-nans solution u ( r , t ) (2:2 response). A different 2-cycle { u (cid:48) , u (cid:48) } , where u (cid:48) = Φ T [ u (cid:48) ] = u (cid:48) ( r ,
0) and u (cid:48) = Φ T [ u (cid:48) ] = u (cid:48) ( r , T ), corresponds to the conduc-tion block solution u (cid:48) ( r , t ) (2:1 response). These threeperiodic solutions are compared in Fig. 3.Stability of time-periodic solutions can be described interms of the dynamics of deviations δ u ( r , t ) = u ( r , t ) − u k ( r , t ) from the reference solution u k . The dynamics of finite disturbances is described by ∂ t δ u = D ∇ δ u + f ( u k + δ u ) − f ( u k ) . (7)For infinitesimal disturbances, (7) reduces to a linearPDE ∂ t δ u = L δ u , (8)where L = D ∇ + J and J ( t ) = ∂ u f | u k ( t ) is the time-dependent Jacobian of the nonlinear term describing cel-lular dynamics evaluated on the periodic orbit.The linear time-evolution operator that advances thesolution of (8) from 0 to t in the tangent space, δ u ( r , t ) = U ( t, δ u ( r ,
0) (9)can be found by integrating (8), which yields a time-ordered exponential U ( t,
0) = exp (cid:26)(cid:90) t L k ( t (cid:48) ) dt (cid:48) (cid:27) . (10)Linear stability of u k ( r , t ) is determined by the eigenval-ues (Floquet multipliers) λ ki of the operator U ( kT, U ( kT, e ki ( r ) = λ ki e ki ( r ) , (11)where e ki ( r ) are the corresponding eigenfunctions (Flo-quet modes). In particular, if | λ ki | < i , anyinfinitesimal disturbance δ u ( r , t ) will eventually decay tozero, so that u ( r , t ) → u k ( r , t ) for t → ∞ . In the follow-ing discussion, we will index the eigenvalues in the orderof decreasing absolute value, | λ k | ≥ | λ k | ≥ | λ k | ≥ · · · .Asymptotic stability, however, does not imply thatinitial disturbances will decay monotonically. If theevolution operator U ( t,
0) is nonnormal (i.e.,
U U † (cid:54) = U † U ), some disturbances may grow transiently, beforeasymptotic decay takes over . As we have shownpreviously in the context of one-dimensional pacedcardiac tissue (Purkinje fibers), such transient growthmay arise, for example, due to the action of closed-loop feedback control, which makes the 1:1 solution lin-early stable. In the present case, the Jacobian J ( t ) isnon-self-adjoint, so both the differential operator L ( t )and the evolution operator U ( t,
0) are nonnormal, hence,generalized linear stability theory is required to de-scribe transient response to perturbations associatedwith changes in the pacing rate.The adjoint of the evolution operator U † ( t,
0) = exp (cid:26)(cid:90) t L † ( t − t (cid:48) ) dt (cid:48) (cid:27) (12)plays an important role in the generalized stability the-ory. It defines the evolution of disturbances δ v in theadjoint of the tangent space, δ v ( r ,
0) = U † ( t, δ v ( r , t ) , (13)backward in time. In practice it is more straightforwardto compute the action of U † ( t,
0) on a vector by solvingthe linear PDE − ∂ t δ v = L † δ v (14)with “initial” condition defined at the final time.Given that time-periodic states u k ( r , t ) could be sta-ble or unstable, we computed them as solutions u k of(6) using a Newton-Krylov solver . The spectrumof the evolution operator U ( t,
0) was computed usingthe Arnoldi method implemented by the Matlab func-tion eigs with accuracy yielding eight significant dig-its. In both instances matrix-free evaluation of the thematrix-vector product U ( t, δ u was accomplished viasimultaneous numerical time-integration of (1) and (8).The matrix-free calculation of the product U ( t, † δ v in-volves backward integration and cannot be accomplishedby integrating (1) and (14) simultaneously. Instead aprecomputed and interpolated solution u k ( r , t ) was usedto compute the solution to (14) using the procedure de-scribed in Ref. 42. The time integrators were imple-mented as OpenCL kernels administered by host codeswritten in C and encapsulated in Matlab mex-files. Thisallowed substantially speeding up the calculations by ex-ecuting them on general purpose Graphics ProcessingUnits (GPUs), at the same time retaining convenient in-terface within Matlab scripts. C. Stability spectra
The 1:1 solution u ( r , t ) exists at all T considered inthis study. It is stable for T > T c ≈ . T < T c . The leading eigenvalue crosses theunit circle along the negative real line, which correspondsto a period doubling bifurcation, at which point the al-ternans solution u ( r , t ) with period 2 T is created. Itcorresponds to discordant alternans with four stationarynodal lines, as Fig. 4(a) illustrates. The alternans iscreated discordant, since the domain is sufficiently large,according to theoretical analysis . In experiments and simulations alternans is typically created con-cordant and becomes discordant as the pacing interval isdecreased. However, this subtle distinction has no bear-ing on the development of conduction block that playsthe crucial role in transition to fibrillation.To further characterize this bifurcation, we computedthe distance (cid:107) u − u (cid:107) between the 1:1 and 2:2 solutionsin the state space, defined using the Euclidean 2-norm (cid:107) u ( r ) (cid:107) = (cid:104) u ( r ) , u ( r ) (cid:105) / , (15)where the scalar product is given by (cid:104) v , u (cid:105) = (cid:90) v ∗ ( r ) · u ( r ) dx dy, (16)and the integral is over the entire computational do-main 0 < x, y < L . The bifurcation is supercritical (a) −15015 (b)
80 84 88 T (ms) k u ! u k T c FIG. 4. The 2:2 solution which corresponds to discordantalternans. (a) The difference (in ms) between two successiveAPDs at T = 80 ms. (b) The distance between the 1:1 and2:2 solutions in the state space. The arrow marks the criticalperiod T c .
60 70 801.271.31.33 T (ms) R unstablestable FIG. 5. Regions in the (
T, R ) parameter plane showing wherethe 2:2 state is stable and where it is not. The solid line isthe neutral stability curve which corresponds to | λ | = 1. (cf. Fig. 4(b)), also in agreement with the theoreticalpredictions , and the 2:2 solution exists for all T < T c (at least down to T = 50 ms).The stability balloon for the 2:2 solution in the ( T, R )plane is shown in Fig. 5. This solution undergoes asubcritical Hopf bifurcation when the neutral stabilitycurve (shown as a solid line) is crossed. The minimum ofthe neutral stability curve is at R c = 1 .
279 and T = 66 . T = 82 ms T = 50 msIm( λ )1 Re( λ )10 (a) Im( λ i )1 Re( λ i )10 (b) Im( λ i )1 Re( λ i )10 (c) FIG. 6. (a) Leading eigenvalue λ shown for 50 ms ≤ T ≤
82 ms in steps of 2 ms. The spectrum { λ i } , i = 1 , , ..., T = 82 ms (b) and T = 50 ms (c). ms. The value R = 1 .
273 considered here is below R c , sothe 2:2 solution turns out to be linearly stable for all T .The corresponding Floquet spectra for different valuesof T are shown in Fig. 6. The change in the leadingeigenvalue λ associated with a decrease in T is shownin Fig. 6(a). Although this eigenvalue never leaves theunit circle, it approaches this circle for T ≈ . R = 1 .
273 and T = 80 ms. In addition, there isone more stable state that exists for the same parametervalues, the one that describes a 2:1 response. For R =1 . T ≈ . T decreases to at least 50ms (that is, over the entire range of T in Fig. 5).Given this information, how can we understand thedifference in the outcomes of the pacing protocols A andB? Both protocols can be viewed as being composedof two stages: the first stage (first 40/50 pacing inter-vals for protocol A/B) determines the initial conditionfor the second stage (pacing with constant T = 80 msthat starts, respectively at t or t ). The perturba-tion δ u B = u B ( r , t ) − u (cf. Fig. 7(c-d)) in proto- (a) 024 (b) −0.200.2 (c) −2024 (d) −0.2−0.100.1 FIG. 7. Comparison of the initial perturbations δ u A ( r ) = u A − u and δ u B ( r ) = u B − u , for the second stage ofprotocols A and B, respectively. (a) and (b) show u - and v -components of δ u A . (c) and (d) show u - and v -componentsof δ u B . col B eventually decays, in agreement with the predic-tion of linear stability analysis, while the perturbation δ u A = u A ( r , t ) − u (cf. Fig. 7(a-b)) in protocolA grows, producing conduction block and wave breakupafter 15 pacing intervals and, eventually, transition topersistent SWC. (Note that in the 2:1 state conductionblock is global, annihilating an entire excitation wave,while in protocol A conduction block is local and leadsto a breakup, but not a total annihilation, of an excita-tion wave.)Since three different stable solutions (2:2, 2:1, andSWC) coexist at T = 80 ms, the asymptotic regime is de-termined by the initial condition. In particular, there is aneighborhood Ω = Ω ∪ Ω of the 2-cycle { u , u } ,such that for all initial conditions u ∈ Ω , the orbit u n approaches this 2-cycle as n → ∞ . Ω is knownas the basin of attraction of the 2-cycle and has a fi-nite size. Similarly, the 2-cycle { u (cid:48) , u (cid:48) } has a basinof attraction Ω = Ω (cid:48) ∪ Ω (cid:48) . In particular, the initialcondition u B ≡ u B ( r , t ) in protocol B lies inside Ω ,so the asymptotic state corresponds to the 2:2 solution.On the other hand, the initial condition u A ≡ u A ( r , t )in protocol A lies outside of both Ω and Ω , so theasymptotic state corresponds to SWC.Since the evolution operator describing the linearizeddynamics is nonnormal, the basin of attraction Ω be-comes quite thin in some directions. Fig. 8 shows atwo-dimensional cartoon of this; the actual dimension-ality of Ω is the same as that of the state space ofthe discretized PDE, 18432 in the present case. We areprimarily interested in the direction in the state space(which corresponds to perturbation shape in the physi-cal space) for which the distance (cid:107) δ u (cid:107) from u to the ! ′Ω q p u u u ! Ω ! Ω ! ′Ω B u A u ! ! δ u optcr FIG. 8. Schematic representation of the basins of attrac-tion of the stable 2T-periodic solutions in a two-dimensionalprojection of the state space. Gray dots represent the initialconditions for the second stage of protocols A and B. Thevector δ u cropt shows the critical linearly optimal perturbation.The white region corresponds to the basin of attraction Ω swc of the SWC state. boundary ∂ Ω of the basin of attraction is the smallest(we will refer to it as the optimal direction), since thedynamics are the most sensitive to perturbations corre-sponding to this direction. The point on ∂ Ω that is theclosest to the 2-cycle defines the critical optimal pertur-bation δ u cropt . Although the shape of Ω , and hence theperturbation δ u cropt , is defined by the nonlinear evolutionoperator (5), as we illustrate below, both the directionof δ u cropt and the strong transient growth associated withperturbations in this direction can be understood reason-ably well within the linear approximation. D. Generalized linear stability analysis
Let us compute how strongly an initial disturbance δ u ( r ,
0) is amplified as a result of transient growth. Withthe help of (9) and (15) we find (cid:107) δ u ( r , t ) (cid:107) = (cid:104) U ( t, δ u ( r , , U ( t, δ u ( r , (cid:105) , = (cid:104) U † ( t, U ( t, δ u ( r , , δ u ( r , (cid:105) . (17)Since U † ( t, U ( t,
0) is self-adjoint, it has real eigenvalues σ i ( t ) ≥ q i ( r , t ), U † ( t, U ( t, q i ( r , t ) = σ i ( t ) q i ( r , t ) , (18)with eigenvalues sorted from largest to smallest. Expand-ing the initial condition in the eigenvector basis yields δ u ( r ,
0) = (cid:88) i η i ( t ) q i ( r , t ) , (19)where η i ( t ) = (cid:104) q i ( r , t ) , δ u ( r , (cid:105) , such that (cid:107) δ u ( r , t ) (cid:107) = (cid:88) i σ i ( t ) η i ( t ) . (20)The maximum of (cid:107) δ u ( r , t ) (cid:107) over all initial disturbanceswith a fixed norm (cid:107) δ u ( r , (cid:107) is reached when all η i ( t ) = 0,except for i = 1, i.e., δ u opt ( r , ∝ q ( r , t ).Let (cid:107) q i ( r , t ) (cid:107) = 1 and define p i ( r , t ) = U ( t, q i ( r , t ).The action of the evolution operator U ( t,
0) can now beexpressed in a simple manner. Multiplying both sides of t=T < FIG. 9. Transient amplification magnitude for T = 80 ms.The singular value σ ( t ) is shown over the discrete set of timeswith step ∆ t = T /
10 = 8 ms. (19) by U ( t,
0) and recognizing that δ u ( r ,
0) is arbitrary,we obtain U ( t,
0) = (cid:88) i p i ( r , t ) σ i ( t ) (cid:104) q i ( r , t ) , ·(cid:105) , (21)which is the singular value decomposition (SVD) of U ( t,
0) with σ i ( t ), p i ( r , t ), and q i ( r , t ) being, respec-tively, the singular values, left singular vectors, and rightsingular vectors of U ( t, q ( r , t ) de-termines the shape of the optimal initial perturbation δ u opt ( r ,
0) that is amplified the most (in the sense ofthe 2-norm (17)) at time t . The left singular vector p ( r , t ) determines the shape this optimal perturbationacquires at time t . Finally, the singular value deter-mines the magnitude of transient amplification σ i ( t ) = (cid:107) δ u opt ( r , t ) (cid:107) / (cid:107) δ u opt ( r , (cid:107) which corresponds to this par-ticular initial condition.For small initial disturbances, the time at which themaximal transient amplification σ ( t ) is achieved is de-termined by the balance between transient growth andasymptotic decay. The weaker the asymptotic decay, thelarger the amplification that can be achieved. Computa-tion of max t σ ( t ) is numerically rather expensive, sinceSVD requires alternate direction integration of the linearPDEs (8) and (14) on the time interval [0 , t ] and thereis no simple relation that allows one to compute the sin-gular values at time t + dt based on their values at time t . We therefore performed this calculation over manypacing periods on a coarse grid with ∆ t = T / σ ( t ) is shown in Fig. 9. Wesee that transient growth is a relatively slow process. Themaximal transient amplification takes place only afteralmost 30 pacing stimuli. At this level of resolution, σ ( t )appears to be a discontinuous function of time. In fact,this is not the case, as the time-resolved calculation overthe interval [26 T, T ] illustrates (cf. Fig. 10). We findthat σ ( t ) has a rather fine structure on a time scale
26 27 28 29 3000.51 t/T σ , k δ u k FIG. 10. Time-resolved transient amplification magnitudeover the interval t ∈ [26 T, T ]. The singular value σ ( t )(solid gray line) and the norm of the perturbation δ u ( r , t ) = U ( t, q ( r , t max ) (dashed black line) for T = 80 ms. Bothquantities were normalized by their maxima. smaller than one pacing period. The maxima of transientamplification have a pattern that is 2 T -periodic; they arelocated at t ≈ (2 n + 0 . T , where n is an integer. Thepeaks’ values, on the other hand, are modulated with aperiod 2 π/ arg( λ ) ≈ . T determined by the leadingFloquet multiplier. The largest transient amplification( σ = 804) is achieved at t max = 28 . T and the secondhighest value is achieved at t = 26 . T . Recall that forprotocol A it takes a time interval of roughly 15 T untilthe dynamics transitions from the neighborhood of the2:2 state to SWC. These times are all of the same orderof magnitude, suggesting that memory effect extends overtimes scales of order 20 T = 1600 ms, much longer thanthe naive estimate (cid:15) − τ u = 250 ms would suggest.The right singular vectors (cf. Fig. 11(a-b)) associatedwith the time instances when σ ( t ) is near the peak value σ ( t max ) are indistinguishable, implying that their spa-tial structure, which determines the optimal initial per-turbation δ u opt ( r ,
0) is a robust feature of the dynamics.Indeed, the time-dependence of σ ( t ) for t (cid:38) T can bereproduced with extremely high accuracy by the norm ofthe disturbance which corresponds to the optimal one,e.g., δ u opt ( r ,
0) = s q ( r , t max ), where s is the amplitudeof the initial disturbance. As Fig. 10 illustrates, the norm (cid:107) δ u opt ( r , t ) (cid:107) is indistinguishable from σ ( t ) once bothhave been normalized by their maximal values. Effec-tively, this means that the initial disturbance s q ( r , t max )is in fact optimal for all times t (cid:38) T . This optimaldisturbance is strongly localized around the pacing site,which means that the dynamics are very sensitive to theperturbation in the timing of the pacing stimuli.The corresponding left singular vector p ( r , t max ) isshown in Fig. 11(c-d)). It is worth noting that, unlikethe right singular vector which is strongly localized, theleft singular vector is not. This spatial delocalization isto a large extent responsible for the apparent large tran-sient amplification with respect to the 2-norm (15), which -3 (a) -8-40 -2 (b) -16-80 × -2 (c) -4-20 -4 (d) -6-303 FIG. 11. Leading singular vectors. (a) and (b) show u -and v -component of the right singular vector q ( r , t max ). (c)and (d) show u - and v -component of the left singular vector p ( r , t max ). Dashed black lines indicate the positions of thenodal lines of the 2:2 solution. involves integration over the entire spatial domain. Fur-thermore, the right singular vector is dominated by the v component, so that the dynamics are more sensitiveto the perturbation of the v variable than the u vari-able. For the left singular vector the opposite is true,hence the perturbation affects the u variable to a muchlarger degree than the v variable. A similar asymmetrywas found for the left and right marginal eigenvectors(response functions and Goldstone modes) of spiral wavesolutions, which similarly describe the cause-and-effectrelationship for small perturbations .As expected, generalized linear stability analysis cor-rectly predicts the dynamics produced by optimal initialdisturbances of sufficiently small amplitude s . Evolvingthe optimal initial perturbation with magnitude s = 10 − using the nonlinear equation (7), we find that δ u opt ( r , t )grows transiently, achieving the predicted amplification σ ( t max ) at time t max (cf. Fig. 12(a)). For t (cid:38) T the asymptotic exponential decay (cid:107) δ u opt ( r , t ) (cid:107) ∝ e κt with the predicted rate κ = ln | λ | / (2 T ) sets in, where | λ | = 0 . s = 10 − leads to the sameresult, with the entire curve simply shifted up, indicatingthat nonlinear effects are still negligible.These results are qualitatively similar to transient am-plification found for traveling wave solutions describingthin fluid films spreading on a solid substrate . In thatcase as well, transient growth of small localized pertur-bations was traced to the nonnormality of the evolu-tion operator . Strong transient amplification of dis-turbances also characterizes traveling excitation wavein models of paced Purkinje fibers with feedback, al- −10 −5 t/T k δ u ( r , t ) k (a) −1 t/T k δ u ( r , t ) k (b) FIG. 12. Transient growth of the optimal initial pertur-bation δ u opt ( r ,
0) = s q ( r , t max ) for different perturbationstrengths s . (a) s = 10 − (black line) and s = 10 − (solidgray line). The dashed line depicts the predicted asymptoticdecay (cid:107) δ u (cid:107) e κt . (b) s = 1 .
227 (blue line), s = 1 .
228 (red line), s = 1 . s = 1 . though in that case the base state corresponds to the 1:1response . In fact, increasing transient amplificationwas identified as the reason for breakdown of feedbackcontrol for long fibers controlled at the pacing site. E. Finite amplitude disturbances
Since the 2:2 state is linearly stable, transition to SWCis subcritical and requires a finite amplitude perturba-tion. The critical magnitude of the perturbation awayfrom u that causes a transition to a different solution( u , u (cid:48) , u (cid:48) , or SWC) is given by the distance from u to the nearest boundary of the corresponding basinof attraction (cf. Fig. 8). For linearly optimal initialperturbations, δ u opt ( r ,
0) = s q ( r , t max ), we find thatthe asymptotic state is u for s ≤ . . ≤ s ≤ . u , which is also a 2:2 state, but with a different phase.For s = 1 . u , butfor s = 1 . T . Increasing s further we findeither persistent SWC or transient SWC that eventually −1 t/T k δ u ( r , t ) k FIG. 13. The norm of the deviation δ u ( r , t ) = u ( r , t ) − u ( r , t ) from the 1:1 solution for initial conditions u ( r ) + s q ( r ; t max ) with s = 1 .
227 (black line) and s = 1 .
228 (grayline). The dashed lines describe the exponential decay (cid:107) δ u (cid:107) ∝ e κ t and exponential growth (cid:107) δ u (cid:107) ∝ e κ t with the rates κ i = ln | λ i | /T predicted by the linearization around u ( r , t ). gives way to one of the 2:1 states ( u (cid:48) or u (cid:48) ).The extreme sensitivity of the outcome to the choiceof initial conditions is a characteristic feature of chaoticattractors that have fractal basin boundaries . Since u and u are parts of the same 2-cycle, both Ω andΩ share a part of their boundary with the basin ofattraction Ω swc of the SWC state. Like a quilt, basinboundaries are composed of a patchwork of pieces, eachone corresponding to the stable manifold of an edge state– a saddle solution embedded into the boundary, whichhas a one-dimensional unstable manifold. The relationbetween optimal perturbations, transient growth, edgestates, and basin boundaries has originally been investi-gated in the context of simple nonlinear PDEs . Morerecently, edge states have been studied extensively asgateways mediating subcritical transition to turbulencein shear fluid flows .The 1:1 state u is an example of such an edge state.It has one unstable direction and its stable manifoldforms the boundary between Ω and Ω . As Fig. 12(b)illustrates, the solutions corresponding to perturbationswith s = 1 .
227 and s = 1 .
228 approach and closely fol-low the same periodic solution for 40 T (cid:46) t (cid:46) T , be-fore eventually separating and approaching, respectively, u and u . This periodic solution is indeed u , asshown in Fig. 13. As expected, the two trajectories ap-proach u ( r , t ) along the least stable direction and sepa-rate along the unstable direction, with the rates predictedby the linearization.The basin boundary ∂ Ω swc of the chaotic attractor liesfurther away from u than the basin boundary ∂ Ω of the 2:2 state u . An upper bound for the distanceto ∂ Ω swc is given by the norm s cr = (cid:107) δ u cropt ( r , (cid:107) =1 . (cid:107) u (cid:107) = O (200) itself. The ratio of the0 (a) (b) (c)(d) (e) (f)(g) (h) (i) FIG. 14. Wave breakup produced by the linearly optimaldisturbance δ u opt ( r ,
0) = s q ( r , t max ) with, s = 1 . u ( r , t ) is shown at (a) t =19 . T , (b) t = 19 . T , (c) t = 19 . T , (d) t = 22 . T , (e) t = 23 . T , (f) t = 23 . T , (g) t = 27 . T , (h) t = 30 . T , and(i) t = 34 T . Dashed white lines indicate the positions of thenodal lines of the 2:2 solution. norms is not as large as the maximal transient amplifica-tion σ ( t max ) = 804, but it is of the same order of mag-nitude, which suggests that linear theory cold be used toestimate the critical disturbance magnitude. However, afully nonlinear calculation is required to obtain a quanti-tatively accurate prediction for both the magnitude andthe shape of the smallest initial disturbance required totrigger a transition from the 2:2 state to SWC.Figure 14 shows the evolution of the initial condi-tion u ( r ,
0) = u + δ u opt ( r , s = 1 . . Eventhough the 2:2 state is linearly stable and the perturba-tion is quite small, transient amplification of the initialdisturbance leads to pronounced alternans which causesconduction block in the lower right corner of the domainat t = t r ≈ . T . Recall that, according to general-ized linear stability theory, the shape of the perturba-tion at the moment when it reaches the largest amplifi-cation is given by the left singular vector of U ( t max , δ u opt ( r , t max ) = σ ( t max ) s p ( r , t max ). The location ofconduction block is indeed found to be consistent withthe shape of the left singular vector p ( r , t max ), which (a) − (b) −0.300.3 FIG. 15. The perturbation δ u opt ( r , t r ) describing conductionblock at t = t r ≈ . T , which corresponds to the initialcondition δ u opt ( r ,
0) = s q ( r , t max ) with s = 1 . u - and v -component, respectively. Dashedblack lines indicate the positions of the nodal lines of the 2:2solution. has maximal amplitude in the same corner of the domain,opposite the pacing site (cf. Fig. 11(c-d)). However,the overall shape of the actual deviation δ u cropt ( r , t r ) = u ( r , t ) − u ( r , t ) (cf. Fig. 15) is noticeably different from p ( r , t max ) (cf. Fig. 11(c,d)).While the first instance of conduction block generatesa pair of spiral waves at t ≈ . T , this does not causean immediate transition to SWC. The spiral waves havetoo little room to sustain themselves: they collide withthe next wavefront emitted by the pacing site and areannihilated at t ≈ T . Another instance of conductionblock occurs at t ≈ . T , again generating two spi-ral waves at t ≈ . T . Now conduction block occurssooner in the cycle and further away from the lower rightcorner of the domain, and the spiral waves survive thecollision with the next wavefront emitted by the pacingsite. (Note, however, that in both instances conductionblock occurs near the nodal line that is the furtherestfrom the pacing site.) Since the period of the spiral wavesis shorter than the pacing period, they gradually invadethe domain, generating further breakups and initiatingsustained SWC.As we discussed previously, from t = t on, pro-tocol A corresponds to constant pacing with T = 80ms. The initial condition u A = u ( r , t ) for the laststage of protocol A corresponds to a fairly significantperturbation away from the 2:2 state u . The norm (cid:107) δ u A (cid:107) = (cid:107) u A − u (cid:107) = 104 . s cr = (cid:107) δ u cropt (cid:107) (andis close to the norm of u itself). On the other hand,the spatial structure of this perturbation (cf. Fig. 7(a-b)) is quite different from that of the optimal pertur-bation (cf. Fig. 11(a-b)) and therefore it does not ex-perience transient amplification, as Fig. 16 illustrates.This difference can be quantified using the scalar prod-uct (cid:104) δ ˆ u A , ˆ q (cid:105) = 0 . δ ˆ u = δ ˆ u / (cid:107) δ ˆ u (cid:107) ), which shows that δ u A hasessentially no component along q . At the same time (cid:104) δ ˆ u A , ˆ p (cid:105) = − . δ u A has a significant com-ponent in the direction of p (cf. Fig. 8).Despite the dramatic difference in the amplitude and1 ( t − t n ) /T k δ u ( r , t ) k FIG. 16. The norm of the perturbations away from the 2:2solution describing the protocols A (black line, n = 40) andB (gray line, n = 50). spatial structure of the perturbations δ u A and δ u cropt ,the sequence of states generated by protocol A (cf. Fig.2(b-d)) looks very similar to that produced by the near-critical optimal perturbation (cf. Fig. 14). One possi-ble explanation is that, before approaching the chaoticattractor, both trajectories closely follow the boundarybetween Ω and Ω swc for an extended period of time(around 20 T for both δ u A and δ u cropt ) and finally theystart to separate from it along the unstable manifold ofthe same edge state that lies on this boundary.Although itself unstable, the basin boundary can befollowed using the bisection algorithm that has beensuccessfully used in the context of fluid turbulence. Forthe system considered here we found that the critical op-timal perturbation δ u cropt approaches a chaotic edge statethat lies on the boundary of Ω . It may in principlebe possible to find nonchaotic edge states embedded into ∂ Ω as well. For instance, the unstable quasiperiodic so-lution born an the subcritical Hopf bifurcation of the 2:2state lies on this boundary and, if it is has one unstableeigenvalue, it would be an example of an edge state. Evensimpler (e.g., time-periodic) edge states may be foundby applying the bisection algorithm to initial conditionsother than u + δ u cropt , however this is far outside thescope of this study.The perturbation that corresponds to protocol B is alsomuch stronger than the critical optimal one, but since itsspatial structure (cf. Fig. 7(c-d)) is also quite differentfrom δ u opt , it does not experience transient amplificationeither, as Fig. 16 illustrates. Since (cid:104) δ ˆ u B , ˆ q (cid:105) = 0 . (cid:104) δ ˆ u B , ˆ p (cid:105) = − . δ u B has essentially no compo-nent along q , but a significant component along p . Thenorm of the perturbation, (cid:107) δ u B (cid:107) = (cid:107) u B − u (cid:107) = 92 . u B to fall withinthe basin of attraction of the 2:2 state (cf. Fig. 8). Hencethe perturbation eventually decays, and the dynamics re-turn to concordant alternans instead of transitioning toSWC. IV. DISCUSSION
The results presented in this paper help improve ourunderstanding of the transition to fibrillation in pacedatrial tissue and the role of alternans. Specifically, theseresults confirm that the first step in the transition causedby an increase in the pacing rate is the development ofdiscordant alternans . Theoretical analysis of alter-nans in one dimension predicts that conduction blockshould occur at the location that is the furthest from thepacing site. On the other hand, experimental and numer-ical evidence suggests that reentry is promoted by steeprepolarization gradients and therefore conductionblock should happen near a nodal line. Indeed, we findthat conduction block that initiates wave breakup andreentry happens to lie close to the nodal line of the sta-ble discordant alternans that is the furtherest from thepacing site. It should be pointed out, however, that thiscorrespondence is not precise, since nodal lines can (anddo) move during the transient evolution from stable dis-cordant alternans to conduction block.However, in a marked departure from the prevailingparadigm, the transition from discordant alternans toconduction block (followed by wave breakup, reentry,and SWC) is not associated with a bifurcation leadingto either an instability or disappearance of the alternansstate. Instead the transition is triggered by strong tran-sient amplification of disturbances away from stable al-ternans. At the pacing interval of 80 ms, certain infinites-imal disturbances were found to be amplified by almostthree orders of magnitude. Transient amplification in-creases as the pacing interval is decreased, reaching amaximum of four orders of magnitude for a pacing in-terval of 66 ms. This is very similar to the trend foundfor linearly stable fluid flows , where transient amplifi-cation increases with the Reynolds number, which playsthe role analogous to the pacing rate.Furthermore, we found that the spatial profile of theoptimal perturbations that experience the strongest tran-sient amplification corresponds to a tiny change in thepacing interval. Hence even small changes in the pacinginterval will generate a strong perturbation away fromthe 2:2 response. For a protocol with a decreasing pac-ing interval, the strength of a perturbation is controlledby the rate at which the pacing interval changes. Sincetransient amplification increases with decreasing pacinginterval, transition from alternans to SWC will happensooner (for a longer pacing interval) if the pacing intervalis decreased faster.Although generalized linear stability analysis is capa-ble of describing transient dynamics qualitatively, it can-not produce a quantitatively accurate description of thetransition, which involves intrinsically nonlinear effects,such as conduction block. Using numerical simulations,we were able to predict that, for a pacing interval of80 ms, transition to SWC is triggered by a linearly op-timal perturbation the norm of which is two orders ofmagnitude smaller than the norm of the alternans state2itself. Even smaller, nonlinearly optimal perturbationsabout the alternans state could trigger a transition toSWC. Computing their magnitude and spatial structure,however, would require a fully nonlinear analysis, fol-lowing, for example, an adjoint-based minimization ap-proach used to characterize transition to turbulence inshear fluid flows .Our results also have interesting implications for theproblem of low-energy defibrillation. Although some ofthe modern approaches have demonstrated that theenergy required for defibrillation could be reduced sub-stantially compared with the standard defibrillation pro-tocols, further dramatic reduction may be possible byexploiting the geometric structure of the state space. Un-derstanding the dynamics in the region of state space sur-rounding an edge state that lies on the basin boundarybetween SWC and time-periodic solutions can help de-sign protocols that can terminate fibrillation with, techni-cally, infinitely small perturbations. In practice, the per-turbation strength will depend on how close the chaoticdynamics approach the basin boundary, but still there isa potential for orders-of-magnitude improvement.The identification of edge states should also help pre-dict the opposite process – when a seemingly tiny pertur-bation to a stable time-periodic rhythm leads to a tran-sition into SWC. More immediately, our results also sup-port an approach designed to prevent fibrillation, ratherthan terminate fibrillation after it started. In that ap-proach feedback control is used to stabilize the normalrhythm, thereby preventing the development of alternansin the first place . This approach should work re-gardless of whether the alternans state is stable or not(e.g., for R > .
279 in the model considered here).The analysis presented here was based on a greatly sim-plified model of cardiac tissue, so while its results mayhave a profound impact on how we view the dynamicalmechanisms that describe the onset of fibrillation, theyneed to be verified using more detailed electrophysiologi-cal models of cardiac tissue in two and three dimensions.In particular, it is important to check whether our mainresults can be reproduced in bi-domain models that arerequired to correctly describe intra- and extra-cellular po-tentials. Finally, it would be interesting to see how thedifference between atria and ventricles (both in termsof the electrophysiology and in terms of dimensionality)affect the dynamical description of transition from dis-cordant alternans to fibrillation.
ACKNOWLEDGMENTS
This material is based upon work supported by theNational Science Foundation under Grant No. CMMI-1028133. The Tesla K20 GPUs used for this researchwere donated by the “NVIDIA Corporation” through theacademic hardware donation program.
REFERENCES D. Mozaffarian, E. J. Benjamin, A. S. Go, D. K. Arnett, et al. ,“Heart disease and stroke statistics-2016 update a report fromthe American Heart Association,” Circulation , e38–e360(2016). S. S. Chugh, R. Havmoeller, K. Narayanan, D. Singh, et al. ,“Worldwide epidemiology of atrial fibrillation: A global burdenof disease 2010 study,” Circulation , 837–847 (2014). J. B. Nolasco and R. W. Dahlen, “A graphic method for the studyof alternation in cardiac action potentials,” J. Appl. Physiol. ,191–196 (1968). D. S. Rosenbaum, L. E. Jackson, J. M. Smith, H. Garam, J. N.Ruskin, and R. J. Cohen, “Electrical alternans and vulnerabilityto ventricular arrythmias,” N. Engl. J. Med. , 235 (1994). D. M. Bloomfield, J. T. Bigger, R. C. Steinman, P. B. Namerow,M. K. Parides, A. B. Curtis, E. S. Kaufman, J. M. Davidenko,T. S. Shinn, and J. M. Fontaine, “Microvolt T-wave alternansand the risk of death or sustained ventricular arrhythmias in pa-tients with left ventricular dysfunction,” Journal of the AmericanCollege of Cardiology , 456–463 (2006). K. Hiromoto, H. Shimizu, Y. Furukawa, T. Kanemori, T. Mine,T. Masuyama, and M. Ohyanagi, “Discordant repolarizationalternans-induced atrial fibrillation is suppressed by verapamil,”Circ J , 1368–1373 (2005). S. M. Narayan, M. R. Franz, P. Clopton, E. J. Pruvot, andD. E. Krummen, “Repolarization alternans reveals vulnerabilityto human atrial fibrillation,” Circulation , 2922–2930 (2011). S. Narayan and J. Zaman, “Mechanistically based mapping ofhuman cardiac fibrillation,” J Physiol , 2399–2415 (2016). R. Majumder, M. C. Engels, A. A. F. de Vries, A. V. Pan-filov, and D. A. Pijnappels, “Islands of spatially discordant APDalternans underlie arrhythmogenesis by promoting electrotonicdyssynchrony in models of fibrotic rat ventricular myocardium,”Scientific Reports , 24334 (2016). K. C. Chang and N. A. Trayanova, “Mechanisms of arrhythmo-genesis related to calcium-driven alternans in a model of humanatrial fibrillation,” Scientific Reports , 36395 (2016). J. M. Pastore, S. D. Girouard, K. R. Laurita, F. G. Akar, andD. S. Rosenbaum, “Mechanism linking T-wave alternans to thegenesis of cardiac fibrillation,” Circulation , 1385–1394 (1999). J. N. Weiss, A. Karma, Y. Shiferaw, P. S. Chen, A. Garfinkel,and Z. Qu, “From pulsus to pulseless: The saga of cardiac alter-nans,” Circulation Research , 1244–1253 (2006). B. Echebarria and A. Karma, “Instability and spatiotemporaldynamics of alternans in paced cardiac tissue,” Phys. Rev. Lett. , 208101 (2002). B. Echebarria and A. Karma, “Amplitude equation approach tospatiotemporal dynamics of cardiac alternans,” Phys. Rev. E ,051911 (2007). P. S. Skardal, A. Karma, and J. G. Restrepo, “Spatiotemporaldynamics of calcium-driven cardiac alternans,” Physical ReviewE , 052707 (2014). F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J. Evans,“Multiple mechanisms of spiral wave breakup in a model of car-diac electrical activity,” Chaos , 852–892 (2002). A. Karma, “Electrical alternans and spiral wave breakup in car-diac tissue,” Chaos , 461–472 (1994). J. J. Fox, R. F. Gilmour Jr, and E. Bodenschatz, “Conductionblock in one-dimensional heart fibers,” Physical Review Letters , 198101 (2002). S. Alonso, M. B¨ar, and B. Echebarria, “Nonlinear physics ofelectrical wave propagation in the heart: a review,” Reports onProgress in Physics , 096601 (2016). E. Surovyatkina, D. Noble, D. Gavaghan, and A. Sher, “Mul-tistability property in cardiac ionic models of mammalian andhuman ventricular cells,” Progress in Biophysics and MolecularBiology , 131–141 (2010). P. S. Skardal and J. G. Restrepo, “Coexisting chaotic and multi- periodic dynamics in a model of cardiac alternans,” Chaos ,043126 (2014). P. Comtois and A. Vinet, “Multistability of reentrant rhythms inan ionic model of a two-dimensional annulus of cardiac tissue,”Physical Review E , 051927 (2005). R. A. Oliver, C. S. Henriquez, and W. Krassowska, “Bistabilityand correlation with arrhythmogenesis in a model of the rightatrium,” Annals of Biomedical Engineering , 577–589 (2005). D. S. Henningson, A. Lundbladh, and A. V. Johansson, “Amechanism for bypass transition from localized disturbances inwall-bounded shear flows,” Journal of Fluid Mechanics , 169–207 (1993). L. N. Trefethen, A. E. Trefethen, S. Reddy, and T. A. Driscoll,“Hydrodynamic stability without eigenvalues,” Science , 578–584 (1993). D. R. Chialvo, D. C. Michaels, and J. Jalife, “Supernormal ex-citability as a mechanism of chaotic dynamics of activation in car-diac purkinje fibers.” Circulation Research , 525–545 (1990). M. B. Rosenbaum, H. H. Blanco, M. V. Elizari, J. O. L´azzari,and J. M. Davidenko, “Electrotonic modulation of the t waveand cardiac memory,” Am. J. Cardiol. , 213–222 (1982). T. J. Hund and Y. Rudy, “Determinants of excitability in cardiacmyocytes: mechanistic investigation of memory effect,” Biophys.J. , 3095–3104 (2000). A. Gizzi, E. M. Cherry, R. F. Gilmour, S. Luther, S. Filippi, andF. H. Fenton, “Effects of pacing site and stimulation history onalternans dynamics and the development of complex spatiotem-poral patterns in cardiac tissue,” Frontiers in Physiology , 71(2013). O. Ziv, E. Morales, Y.-k. Song, X. Peng, K. E. Odening, A. E.Buxton, A. Karma, G. Koren, and B.-R. Choi, “Origin of com-plex behaviour of spatially discordant alternans in a transgenicrabbit model of type 2 long QT syndrome.” The Journal of phys-iology , 4661–4680 (2009). G. Byrne, C. D. Marcotte, and R. O. Grigoriev, “Exact coherentstructures and chaotic dynamics in a model of cardiac tissue,”Chaos , 033108 (2015). C. D. Marcotte and R. O. Grigoriev, “Unstable spiral waves andlocal Euclidean symmetry in a model of cardiac tissue,” Chaos , 063116 (2015). F. H. Fenton and E. M. Cherry, “Models of cardiac cell,” Schol-arpedia , 1868 (2008). Z. Qu, A. Garfinkel, P. S. Chen, and J. N. Weiss, “Mechanismsof discordant alternans and induction of reentry in simulatedcardiac tissue,” Circulation , 1664–1670 (2000). S. A. Orszag and A. T. Patera, “Subcritical transition to turbu-lence in plane channel flows,” Phy. Rev. Lett. , 989 (1980). H. Chat´e and P. Manneville, “Transition to turbulence via spatio-temporal intermittency,” Phys. Rev. Lett. , 112 (1987). A. Garz´on, R. O. Grigoriev, and F. H. Fenton, “Model-basedcontrol of cardiac alternans in Purkinje fibers,” Phys. Rev. E ,041927 (2011). A. Garz´on, R. O. Grigoriev, and F. H. Fenton, “Continuous-timecontrol of alternans in long purkinje fibers,” Chaos , 033124(2014). P. J. Ioannou and B. F. Farrell, “Generalized stability theorypart 1: Autonomous operators,” J. Atmos. Sci. , 2025–2040(1996). P. J. Ioannou and B. F. Farrell, “Generalized stability theorypart 2: Nonautonomous operators,” J. Atmos. Sci. , 2041–2053 (1996). R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly re-started Arnoldi iteration,” SIAM J. Matr. Anal.Appl. , 789 (1996). C. D. Marcotte and R. O. Grigoriev, “Adjoint eigenfunctions oftemporally recurrent single-spiral solutions in a simple model ofatrial fibrillation,” Chaos , 093107 (2016). C. Marcotte and R. O. Grigoriev, “Implementation of PDEmodels of cardiac dynamics on GPUs using OpenCL,” (2012), arXiv:1309.1720 . M. L. Walker and D. S. Rosenbaum, “Repolarization alternans:implications for the mechanism and prevention of sudden cardiacdeath,” Cardiovasc. Res. , 599 (2003). M. A. Watanabe, F. H. Fenton, S. J. Evans, H. M. Hastings, andA. Karma, “Mechanisms for discordant alternans,” J. Cardiovasc.Electrophysiol. , 196–206 (2001). A. L. Bertozzi and M. P. Brenner, “Linear stability and transientgrowth in driven contact lines,” Physics of Fluids , 530–539(1997). R. O. Grigoriev, “Transient growth in driven contact lines,”Physica D , 105–116 (2005). S. W. McDonald, C. Grebogi, E. Ott, and J. A. Yorke, “Fractalbasin boundaries,” Physica D , 125–153 (1985). A. Handel and R. O. Grigoriev, “Transient dynamics and non-linear stability of spatially extended systems,” Phys. Rev. E ,036302 (2006). T. M. Schneider, B. Eckhardt, and J. A. Yorke, “Turbulencetransition and the edge of chaos in pipe flow,” Phys. Rev. Lett. , 034502 (2007). J. F. Gibson, J. Halcrow, and P. Cvitanovi´c, “Visualizing thegeometry of state-space in plane Couette flow,” J. Fluid Mech. , 107–130 (2008). T. Konta, K. Ikeda, M. Yamaki, K. Nakamura, K. Honma,I. Kubota, and S. Yasui, “Significance of discordant st alternansin ventricular fibrillation.” Circulation , 2185–2189 (1990). H. Tachibana, I. Kubota, M. Yamaki, T. Watanabe, and H. To-moike, “Discordant s-t alternans contributes to formation of reen-try: a possible mechanism of reperfusion arrhythmia,” Am. J.Physiol. , H116 (1998). A. Meseguer and L. N. Trefethen, “Linearized pipe flow toreynolds number 10 ,” Journal of Computational Physics ,178–197 (2003). C. C. T. Pringle and R. R. Kerswell, “Using nonlinear transientgrowth to construct the minimal seed for shear flow turbulence,”Phys. Rev. Lett. , 154502 (2010). R. A. Winkle, R. H. Mead, M. A. Ruder, V. Gaudiani, W. S.Buch, B. Pless, M. Sweeney, and P. Schmidt, “Improved lowenergy defibrillation efficacy in man with the use of a biphasictruncated exponential waveform,” American Heart Journal ,122–127 (1989). F. H. Fenton, S. Luther, E. M. Cherry, N. F. Otani, V. Krinksy,A. Pumir, E. Bodenschatz, and R. F. G. Jr., “Termination ofatrial fibrillation using pulsed low-energy far-field stimulation,”Circulation , 467–476 (2009). S. Luther, F. H. Fenton, B. G. Kornreich, A. Squires, P. Bittihn,D. Hornung, M. Zabel, J. Flanders, A. Gladuli, L. Campoy, E. M.Cherry, G. Luther, G. Hasenfuss, V. I. Krinsky, A. Pumir, R. F.Gilmour, and E. Bodenschatz, “Low-energy control of electricalturbulence in the heart,” Nature , 235–9 (2011). B. Echebarria and A. Karma, “Spatiotemporal control of cardiacalternans,” Chaos , 923–930 (2002). S. Dubljevic, S.-F. Lin, and P. D. Christofides, “Studies of feed-back control of cardiac alternans,” Comp. Chem. Eng.32