Unidirectional Pinning and Hysteresis of Spatially Discordant Alternans in Cardiac Tissue
aa r X i v : . [ n li n . PS ] J a n Unidirectional Pinning and Hysteresis of Spatially Discordant Alternans in Cardiac Tissue
Per Sebastian Skardal, ∗ Alain Karma, and Juan G. Restrepo Department of Applied Mathematics, University of Colorado at Boulder, Colorado 80309, USA Physics Department and Center for Interdisciplinary Research on Complex Systems, Northeastern University, Boston, MA 02115, USA
Spatially discordant alternans is a widely observed pattern of voltage and calcium signals in cardiac tissuethat can precipitate lethal cardiac arrhythmia. Using spatially coupled iterative maps of the beat-to-beat dynam-ics, we explore this pattern’s dynamics in the regime of a calcium-dominated period-doubling instability at thesingle cell level. We find a novel nonlinear bifurcation associated with the formation of a discontinuous jumpin the amplitude of calcium alternans at nodes separating discordant regions. We show that this jump unidi-rectionally pins nodes by preventing their motion away from the pacing site following a pacing rate decrease,but permitting motion towards this site following a rate increase. This unidirectional pinning leads to stronglyhistory-dependent node motion that is strongly arrhythmogenic.
PACS numbers: 87.19.Hh, 05.45.-a,89.75.-k
The study of period-two dynamics in cardiac tissue hasbecome an important topic of research in the physics [1]and biomedical communities [2]. The term alternans de-scribes beat-to-beat alternations of both action potential du-ration (APD) and peak intracellular calcium concentration( Ca peaki ). Heart cells generically exhibit alternans when theyare paced rapidly or in pathological conditions. Interest inalternans during the last decade has stemmed from the dis-covery that APD alternations can become “spatially discor-dant” in tissue [3], meaning that APD alternates with oppo-site phases in different regions [4, 5]. Spatially discordant al-ternans (SDA) dynamically creates spatiotemporal dispersionof the refractory period during which cells are not excitable,thereby promoting wave blocks and the onset of lethal cardiacarrhythmias [2].To date, our theoretical understanding of SDA is well devel-oped for the case where APD alternans results from an insta-bility of membrane voltage ( V m ) dynamics at the single celllevel, which originates from the restitution relation betweenAPD and the preceding diastolic interval (DI) between twoaction potentials. Numerical simulations [6] have shown that“nodes”, which are line defects with period-1 dynamics sep-arating discordant regions of period-2 oscillations of oppositephases, can form spontaneously in paced homogeneous tissuedue to conduction velocity (CV) restitution, the relationshipbetween action potential propagation speed cv and DI. In ad-dition, node formation has been understood theoretically inan amplitude equation framework [7, 8] to result from a pat-tern forming linear instability that amplifies spatially periodicstationary or traveling modulations of alternans amplitude.Despite this progress, our theoretical understanding of SDAremains incomplete. Both experiments [9] and ionic modelsimulations [10] have shown that Ca peaki can alternate evenwhen V m is forced to be periodic with a clamped action po-tential waveform, demonstrating that alternans can also resultfrom an instability of intracellular calcium dynamics. Al-though alternans are presently believed to be predominantly Ca i -driven in many instances, our understanding of nodes inthis important case remains limited. Numerical simulationshave shown a qualitatively similar role of CV-restitution in SDA formation for V m - and Ca i -driven alternans [4, 11], butmore complex behaviors for the latter case depending on thestrength of Ca i -driven instability [12] and the nature of Ca i - V m coupling [11, 13]. However, a theoretical framework to in-terpret both computational and experimental observations hasremained lacking.In this Letter, we extend the theoretical framework of [7] touncover novel aspects of SDA formation for Ca i -dominatedinstability and validate our theoretical predictions with de-tailed ionic model simulations. A major finding is that nodemotion can be pinned in one direction owing to the formationof a discontinuous jump in calcium alternans amplitude at anode where only V m exhibits period-1 dynamics. This jumpleads to strongly history-dependent SDA evolution and alsoalters fundamentally the spacing between nodes. We summa-rize here our main results and additional details of both theoryand simulations will be described elsewhere.We start our analysis from the system of spatially coupledmaps of the general form C n +1 ( x ) = f c [ C n ( x ) , D n ( x )] , (1) A n +1 ( x ) = Z L G ( x, x ′ ) f a [ C n +1 ( x ′ ) , D n ( x ′ )] dx ′ , (2)where A n ( x ) , D n ( x ) , and C n ( x ) denote the APD, DI, and Ca peaki , respectively, at beat n and position x , and G ( x, x ′ ) captures the cumulative effect of electrotonic ( V m -diffusive)coupling during one beat. For a cable of length L paced at x = 0 with no flux boundary conditions on V m at both ends, G ( x, x ′ ) = G ( x − x ′ ) + G ( x + x ′ ) + G (2 L − x − x ′ ) ,with G ( x ) = H ξ ( x ) h wxξ (cid:16) − x ξ (cid:17)i where H ξ is Gaus-sian with standard deviation ξ (see Appendix B of [7]) and ξ = √ D V APD ∗ and w = 2 D V /cv ∗ are two intrinsic length-scales expressed in terms of the APD and CV at the alter-nans bifurcation (APD ∗ and cv ∗ , respectively), and D V isthe diffusion constant of V m in the standard cable equation ˙ V m = D V ∂ x V m − I ion . Furthermore, CV-restitution causesthe activation interval T n ( x ) ≡ A n ( x ) + D n ( x ) to vary frombeat to beat along the cable as [7, 14] T n ( x ) = τ + Z x dx ′ cv ( D n ( x ′ )) − Z x dx ′ cv ( D n − ( x ′ )) , (3)where τ is the imposed period at the paced end ( x = 0 ).To complete the model, we need to specify the forms of f a and f c . Since we are interested in understanding the genericbehavior of alternans, we choose simple phenomenologicalforms of those maps defined implicitly by f c /C ∗ = 1 − rc n + c n + αd n , (4) f a /A ∗ = 1 + βd n + γc n +1 , (5)where c n ≡ ( C n − C ∗ ) /C ∗ , d n = ( D n − D ∗ ) /A ∗ , and wealso define a n = ( A n − A ∗ ) /A ∗ . With this choice C n = C ∗ , A n = A ∗ , and D n = D ∗ are trivial fixed points correspond-ing to c n = a n = d n = 0 . Moreover, c n , a n , and d n mea-sure the departure of Ca peaki , APD, and DI from those fixedpoint values during alternans. The cubic polynomial in c n inEq. (4) models a period-doubling bifurcation of the intracellu-lar calcium dynamics with the amplitude of Ca peaki alternansincreasing with the degree of calcium instability r . The term d n in Eq. (5) incorporates APD-restitution. The other crossterms in Eqs. (4) and (5) model the bi-directional V m − Ca i coupling taken to be positive in both directions ( α > and γ > ), corresponding to the typical case of locally in-phaseAPD and Ca peaki alternans.To complete the derivation of maps describing the dynam-ics of a n , d n , and c n , we linearize Eq. (3) about the fixedpoint, which yields a n + d n = − Z x dx ′ Λ ( d n − d n − ) / (6)where Λ ≡ cv ∗ / (2 cv ′∗ ) . This linearization is formally validas long as the amplitude of DI alternans induced by Ca i alter-nans is small enough that we can locally neglect the curvatureof the CV-restitution curve, or | ( r − / γA ∗ cv ′′∗ /cv ′∗ | ≪ .Furthermore, we assume that the evolution of alternans ampli-tude is sufficiently slow that we can make the approximation d n − ≈ − d n . This assumption is valid close to the alternansbifurcation. Substituting d n − = − d n in Eq. (6) and differen-tiating both sides, we obtain a differential equation for d n ( x ) that can be solved exactly. Because the pacing rate is fixed at x = 0 , giving a n (0) + d n (0) = 0 , this yields d n ( x ) = − a n ( x ) + e − x/ Λ Z x dx ′ Λ e x ′ / Λ a n ( x ′ ) . (7)The dynamics is completely specified by Eq. (7) together withthe maps obtained by inserting Eqs. (4)-(5) into Eqs. (1)-(2): c n +1 ( x ) = − rc n ( x ) + c n ( x ) + αd n ( x ) , (8) a n +1 ( x ) = Z L G ( x, x ′ ) [ βd n ( x ′ ) + γc n +1 ( x ′ )] dx ′ . (9) Λ r ( Λ )r ( Λ )(a) c > 0c = 0 (b) (c) smooth discontinuousc > 0 c(x)a(x)(b) 0 5 10 15 20−1−0.500.51 x c(x)a(x)(c)20 40 60 80 100−1−0.500.51 cell n c(n)a(n)(d) 20 40 60 80 100−1−0.500.51 cell n c(n)a(n)(e) FIG. 1: (Color online) (a) Nature of steady-state solutions in the r - Λ plane. From left to right, no alternans ( c = 0 ), smooth calcium pro-files ( c > , smooth), and discontinuous calcium profiles ( c > , dis-continuous). (b) Smooth traveling and (c) discontinuous stationary c ( x ) (blue) and a ( x ) (dashed red) profiles from simulating Eqs. (7)-(9), using ( r, Λ) = (0 . , and (1 . , , respectively. Alternansprofiles obtained from a detailed ionic model [15] analogously show-ing (d) smooth traveling and (e) discontinuous stationary solutions. Note that the pacing rate τ no longer appears in the final equa-tions but is still contained implicitly in the fact that D ∗ , andhence the CV-restitution slope and Λ , can depend on τ . Also,since d n − = − d n in steady-state, Eqs. (7)-(9) remain validin steady-state even further from the bifurcation.In Fig. 1, we present the results of different alternans behav-ior obtained from a numerical survey of Eqs. (7)-(9) where wevary systematically CV-restitution, which becomes shallowerwith increasing Λ , and the strength of Ca i -driven instabilitythat increases with r . In Fig. 1(a) we summarize the natureof steady-state solutions in this parameter space. Small r val-ues yield no alternans solutions ( c = 0 ) where both steady-state c ( x ) and a ( x ) are identically zero. When r is increasedwe find a first bifurcation at a value r (Λ) where steady-statesolutions for both c ( x ) and a ( x ) become non-zero and formsmooth waves ( c > , smooth). If the asymmetry of G , givenby w , is not too small these waves are stationary, otherwisethey move towards the pacing site with a constant velocity,as in the voltage dominated case [7]. For all work presentedhere w = 0 was used, yielding traveling waves in the smoothregime. We found qualitatively similar results for positive w .When r is increased further we find a second bifurcation ata value r (Λ) where calcium alternans profiles become sta-tionary and discontinuous at the nodes separating out of phaseregions ( c > , discontinuous) while a ( x ) remains smoothdue to the smoothing effect of voltage diffusion. Exampleprofiles of steady-state c ( x ) (blue dots) and a ( x ) (dashed red)are shown in Figs. 1(b) and (c) from the smooth and discon-tinuous regions, using ( r, Λ) = (0 . , and (1 . , , re-spectively. For all figures presented in this Letter we use pa-rameters α = γ = √ . , β = 0 , ξ = 1 , and w = 0 . Forcomparison, in Figs. 1(d) and (e) we show c ( n ) and a ( n ) pro-files inferred from numerical simulation of the detailed ionicmodel in Ref. [15], where n indexes individual cells, usingparameter values that give smooth traveling profiles and dis-continuous stationary profiles, respectively. Traveling profileshave arrows indicating movement.The onset of alternans at r (Λ) is mediated by an abso-lute instability analogous to that studied in Ref. [7] for thevoltage-driven case. For β = 0 a linear stability analysisyields thresholds of r (Λ) = 1 − η + 3 ηξ / / (4Λ / ) and − η + ξ ( w Λ) − for the instability of the traveling and sta-tionary modes, respectively, where η = αγ . Furthermore, thewavelength at onset is πξ / Λ / / √ and π ( w Λ) / in thetraveling and stationary cases, respectively, which agrees withthe voltage-driven case in Ref. [7]. Similar expressions can beobtained for β = 0 . Numerical simulations (not shown) are ingood agreement with these theoretical results.We now concentrate on the discontinuous regime that isthe primary focus of this letter. To characterize calcium al-ternans profiles in this regime [cf. Fig. 1(c)], we examinefirst stationary steady-state period-two profiles and substitute c ( x ) = c n ( x ) = − c n +1 ( x ) into Eq. (8). After differentiatingEq. (8) with respect to x and some manipulations, we obtain Λ c ′ ( x ) = c ( x ) − ( r − c ( x ) − α Λ a ′ ( x )( r − − c ( x ) . (10)Thus, when alternans grow from c ∼ with r > r (Λ) ,if c ( x ) = c − = ± p ( r − / the derivative diverges and c ( x ) becomes discontinuous. Through the discontinuity, thequantity αd ( x ) in Eq. (8) remains smooth, so finding theother root of the cubic ( r − c = c + αd gives the valueof c ( x ) at the latter end of the discontinuity. This gives c ( x ) = c + = ∓ p ( r − / and a total jump of amplitude | c + − c − | = p r − . To measure the asymmetry at a node,we introduce the quantity ∆ ≡ || c + | − | c − || / p ( r − / .We will refer to a discontinuity where c ( x ) jumps from c − = ± p ( r − / to c + = ∓ p ( r − / as a normal jump . Aremarkable property of this jump is that the limiting values c + and c − on either side of the node depend only on the strength r of Ca i -driven instability, and is independent of all the otherparameters Λ , η , β , ξ , and w . Experimentally, r = 1 is thepoint in parameter space where an isolated myocyte pacedwith a periodic AP-clamp waveform, or a tissue paced at onepoint with negligible CV-restitution ( Λ = ∞ ) bifurcates to al-ternans. Hence, the ratio c + /c − can be used to deduce r intissue experiments or simulations under a finite effect of CV-restitution, and hence to relate single-cell and tissue behavior. Λ | c + / − | |c − ||c + |00.51 ∆ FIG. 2: (Color online) Using r = 1 . , | c − | (solid blue) and | c + | (dashed red) versus Λ from an initial profile with Λ = 10 . Inset: ∆ . When starting from the unstable base solution without al-ternans ( a = c = 0 ) in the regime r > r (Λ) , SDA formsdynamically as a periodic pattern of discontinuous nodes withnormal jumps. A unique feature of SDA evolution in thisregime, which is entirely absent for V m -dominated instabil-ity, is that both the node positions and alternans profiles candepend strongly on the history of how the parameters r and Λ are varied. If Λ or r are increased starting from a profilewith normal jumps, the position of the nodes remains con-stant, but the shape of the profile deforms in such a way thatthe jump in Ca i -alternans profile becomes symmetrical aboutthe node, i.e. both | c − | and | c + | approach the same limitingvalue | c ± | = √ r − where ∆ vanishes. This shows that,if initial conditions contain discontinuous nodes, jumps neednot be normal in steady-state if c ( x ) = c − = ± p ( r − / isnot attained. This is shown in Fig. 2 where we plot | c − | (solidblue) and | c + | (dashed red) using r = 1 . as we increase Λ from . As Λ is increased, | c − | and | c + | tend toward one an-other. We superimpose theoretical values of | c ± | for the nor-mal jump and Λ − = 0 cases in dot-dashed black for compari-son, noting that | c − | , | c + | vary smoothly between these valuesfor intermediate values of Λ . In the inset we plot ∆ versus Λ ,noting that ∆ → as Λ → ∞ . If Λ is then decreased back un-til it reaches its original value (not shown), the profile recoversits original shape. However, if Λ or r are decreased startingat a point where the jumps are normal, the pattern close tothe node preserves its shape, but the node moves towards thepacing site. Importantly, if Λ or r are increased back after thenode has moved, the node does not return to its original posi-tion, but rather its shape will deform to become symmetricalas described above. Since no parameter change can inducethe node to move away from the pacing site, node motion is unidirectionally pinned . We note that we have also observedunidirectional pinning in our ionic model simulations [15].When the node is unpinned, we find that the location ofthe first node, denoted x , scales linearly with Λ , suggestingthat the node spacing is independent of electrotonic coupling.This linear scaling with Λ in the discontinuous regime is tobe contrasted with the scaling of the node spacing for smoothalternans profiles (e.g. ξ / Λ / for w = 0 ), which dependsstrongly on electrotonic coupling. Physically, this linear scal-ing reflects the fact that electrotonic coupling has a negligi-ble effect on the outer scale where the alternans profile variesslowly on a scale ∼ Λ , and only becomes relevant on a scale ∼ ξ near the nodes. This only adds a subdominant correctionof order ξ to the x ∼ Λ scaling. Mathematically, it can be Λ (a) (i) (i)(ii) (ii) 1 2 3−101 x c ( x ) (b) (i)(initial)(ii)0.20.40.60.81 ∆ ∆ x (c) path (i) 0.20.40.60.81 x ∆ x (d) path (ii) FIG. 3: (Color online) (a) Paths (i) and (ii) (solid and dashed) in ( r, Λ) . (b) Initial profile (solid black) and final profiles (solid anddashed blue) after moving along path (i) and (ii). (c), (d) Asymmetry ∆ (circles) and x (crosses) along paths (i) and (ii), respectively. related to the fact that Λ scales out of Eq. (10) in the limit Λ ≫ ξ if one uses the scaled variable ˜ x = x/ Λ instead of x .To investigate the consequences of unidirectional pinning,we investigated the pattern evolution in response to multi-ple parameter changes. Namely, we changed r and Λ fol-lowing two different paths that connect the same points in ( r, Λ) . Starting with a pattern with normal jumps at ( r, Λ) =(1 . , , we move to (1 . , first by increasing r andthen decreasing Λ [path (i)] and vice versa [path (ii)]. De-spite the same start and end parameters, the resulting profilecharacteristics vary significantly depending on the path fol-lowed, as shown in Fig. 3. In Fig. 3(a) paths (i) and (ii) aredenoted by solid and dashed blue arrows, respectively, withthe start and endpoints denoted as a black circle and square.In Fig. 3(b) we zoom in on the first node of the initial pro-file (solid black curve) and final profiles after moving alongpath (i) (solid blue) and (ii) (dashed blue). Consistent withour summary above, x remains constant along path (i) while ∆ decreases as r is increased, after which ∆ increases as Λ isdecreased. Along path (ii) x decreases with Λ , then remainsconstant while ∆ decreases as r is increased. This is shown inFig. 3 (c) and (d), which show ∆ (black circles) and x (redcrosses) measured along paths (i) and (ii), respectively.In conclusion, we have extended our basic theoreticalunderstanding of SDA dynamics to the important case ofcalcium-driven instability. Furthermore we have made a num-ber of new experimentally testable predictions, which we havevalidated by detailed ionic model simulations. The main pre-diction is that node motion becomes unidirectionally pinnedwhen the Ca i alternans profile becomes spatially discontinu-ous above a threshold of Ca i -driven instability. This predic-tion could be tested by first increasing progressively the pac-ing rate (decreasing the inverse restitution slope Λ ), therebycausing the node to move towards the pacing site, as has al-ready been observed in some experiments on APD SDA [4, 5],and then decreasing the pacing rate to its original value. If the Ca i -alternans profile exhibits a jump at the node, thenode should remain stationary. Increasing pacing rate canalso cause several parameters to change, including the de-gree of Ca i -driven instability. However, we have shown thathistory-dependent SDA evolution is robust to multiple param-eter changes (Fig. 3), and hence should be observable in morecomplex situations. We emphasize that unidirectional pin-ning is a purely dynamical phenomenon independent of in-trinsic tissue heterogeneities, which can also potentially pinnode motion. However, we expect pinning due to tissue het-erogeneities to be generally bi-directional, and hence distin-guishable from unidirectional dynamically-induced pining. Asecond prediction is that the spatial jump in Ca i -alternans am-plitude displays remarkably universal features. The magni-tude and asymmetry of this jump are insensitive to most pa-rameters except the degree of Ca i -driven instability, and bothquantities are generally history-dependent.Unidirectional pinning generally makes it harder to elimi-nate SDA by node motion once they are formed. We there-fore expect SDA to be more arrhythmogenic for Ca i - than V m -dominated instability. Given that alternans are believedto be predominantly Ca i -driven in common pathologies suchas heart failure, SDA may play an even more important rolethan previously thought in such pathologies.The work of A.K. and J.G.R. was supported in part by Na-tional Institutes of Health grant No. P01 HL078931. ∗ Electronic address: [email protected][1] A. Karma and R. F. Gilmour, Physics Today , 51 (2007).[2] J. N. Weiss et al., Circ. Res , 1244 (2006); Weiss et al., Circ.Res. , 98 (2011).[3] J. M. Pastore et al., Circulation, , 1385 (1999).[4] H. Hayashi et al., Biophys. J. (2), 448 (2007).[5] O. Ziv et al., J. Physiol. (19), 4661 (2009); S. Mironov, J.Jalife, and E. G. Tolkacheva, Circulation, , 17 (2008).[6] M. A. Watanbe et al., J. Cardiovasc. Electrophysiol. (2), 207(2001); Z. Qu et al., Circulation , 1664 (2000).[7] B. Echebarria and A. Karma, Phys. Rev. Lett. , 208101(2002); Phys. Rev. E , 051911 (2007).[8] D. G. Schaeffer and S. Dai, SIAM J. Appl. Math , 704 (2008).[9] E. Chudin et al., Biophys. J. , 2930 (1999).[10] Shiferaw et al., Biophys J. , 3666 (2003); J. G. Restrepo, J.N. Weiss, and A. Karma, Biophys. J. (8), 3767 (2008).[11] D. Sato et al., Circ. Res. , 520 (2006).[12] D. Sato et al., Biophys. J. , 33 (2007).[13] X. Zhao, Phys. Rev. E , 011902 (2008).[14] M. Courtemanche, L. Glass and J. P. Keener, Phys. Rev. Lett. , 2182 (1993).[15] A. Mahajan et al., Biophys. J. , 392 (2008). The parametersthat have qualitatively similar effects as r and Λ are the SRrelease slope u and a parameter α that changes the rate of re-covery from inactivation of the sodium current I Na ( τ j → ατ j where τ j is the time scale for the j -gate), respectively. Figs. 1(d)and (e) used u = 22 . and . ms − , respectively, α = 1 ,BCL = 340 ms, D V / ∆ x = 0 . ms −1