Memory effects, transient growth, and wave breakup in a model of paced atrium
The mechanisms underlying cardiac fibrillation have been investigated for over a century, but we are still finding surprising results that change our view of this phenomenon. The present study focuses on the transition from normal rhythm to atrial fibrillation associated with a gradual increase in the pacing rate. While some of 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 model of 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 than a conventional instability. It is mathematically very similar to subcritical, or bypass, transition from laminar fluid flow to turbulence, which allows many of the tools developed in the context of fluid turbulence to be used for improving our understanding of cardiac arrhythmias.
Atrial fibrillation is a common type of cardiac arrhythmia characterized by spatially uncoordinated high-frequency patterns of electrical excitation waves. Several different explanations of how fibrillation arises from a highly coordinated and regular normal rhythm have been offered. The most common is a dynamical scenario that involves several stages. First, a milder type of arrhythmia known as alternans, characterized by regular spatial and temporal modulation in the duration of the excitation, is created as a result of a period doubling or Hopf bifurcation. In the second stage, the modulation grows until some portion of an excitation wave fails to propagate (which is known as conduction block), producing wave breakup. In later stages, spiral waves segments, or wavelets, form and undergo subsequent breakups until a complicated dynamical equilibrium is established, with wavelets continually annihilating and reemerging. The current paradigm is that each of these stages represents a separate bifurcation that leads to either a change in the stability, or the disappearance, of the solution describing a previous stage. Our analysis of a simple model of paced cardiac tissue shows that transitions between different stages can happen without any bifurcations taking place. In particular, conduction block can be caused by strong transient amplification of finite amplitude disturbances, which is a flip side of a very well-known memory effect that describes the sensitivity of the asymptotic dynamics of paced cardiac tissue to the entire pacing history.
Fibrillation, a state of uncoordinated contraction of the heart muscle, can occur in the ventricles, leading to sudden cardiac arrest, a major cause of death Mozaffarian et al. (2016), or the atria, where it creates significant health complications Chugh et al. (2014). Understanding of the transition from the normal cardiac rhythm, also known as the 1:1 response, to fibrillation is critical for the development of therapies for the prevention or termination of this dangerous arrhythmia. Abundant experimental evidence points to the role of alternans Nolasco and Dahlen (1968), a beat-to-beat alternation in the action potential duration (APD) at high excitation rates, also known as the 2:2 response, as a precursor of fibrillation both in the ventricles Rosenbaum et al. (1994); Bloomfield et al. (2006) and in the atria Hiromoto et al. (2005); Narayan et al. (2011); Narayan and Zaman (2016). Strong correlation between alternans and fibrillation was also found in the simulations employing high-fidelity models of two-dimensional ventricular Majumder et al. (2016) and atrial Chang and Trayanova (2016) tissue.
In tissue, alternans manifests as spatial modulation in the width of the excitation waves Pastore et al. (1999), which leads to dynamical variation in the tissue refractoriness. In particular, spatially discordant alternans, in which APD alternates out of phase in different regions of the heart, markedly enhances dispersion of refractoriness, increasing the likelihood of conduction block and, in two- or three-dimensional tissue, reentry Weiss et al. (2006). Theoretical analysis of propagating excitation waves in one dimension based on amplitude equations Echebarria and Karma (2002a, 2007) show that voltage-driven alternans are described by stable periodic or quasiperiodic solutions, with the tissue size and dispersion of conduction velocity determining whether the alternans is concordant or discordant. The same is true of calcium-driven alternans Skardal, Karma, and Restrepo (2014), which qualitatively behave like their voltage-driven cousins. Theoretical studies of alternans in higher dimensions Echebarria and Karma (2007) are too limited in scope to make any specific predictions regarding the transition from alternans to fibrillation, but based on experiments and numerical simulations Fenton et al. (2002) we know that conduction block that leads to wave breakup plays a key role.
As the pacing rate or the tissue size in increased, conduction block is predicted to arise when alternans becomes unstable leading to unbounded growth of the amplitude of modulation Karma (1994); Echebarria and Karma (2002a); Fox, Gilmour Jr, and Bodenschatz (2002). Conduction block leads to a time-periodic response such as 2:1 (every other excitation wave is blocked) in one dimension, but typically generates wave breakup and transition to spiral wave chaos (SWC) in two dimensions or scroll wave chaos in three dimensions (respectively, atrial and ventricular fibrillation) Alonso, Bär, and Echebarria (2016). As we show in the present study, there is another possibility: conduction block can be triggered by very small, but finite, perturbations of stable discordant alternans.
Although this result appears counterintuitive at first, it is not completely unexpected. It is well-known that paced cardiac tissue displays multistability in zero Surovyatkina et al. (2010), one Skardal and Restrepo (2014), and two Comtois and Vinet (2005); Oliver, Henriquez, and Krassowska (2005) dimensions. Transitions between different stable attractors in spatially extended dissipative nonequilibrium systems are in fact not uncommon. Transition from stable laminar flow to turbulence in shear fluid flows is perhaps the oldest and best known example. In the case of fluid turbulence, this is known as a bypass transition Henningson, Lundbladh, and Johansson (1993) and involves strong transient amplification of small disturbances Trefethen et al. (1993).
Although no direct evidence of the role of transient amplification in the dynamics of cardiac tissue is available at present, indirect evidence is provided by the so-called memory effect Chialvo, Michaels, and Jalife (1990); Rosenbaum et al. (1982); Hund and Rudy (2000). This term refers to the dependence of asymptotic states of paced tissue preparations on the entire pacing history rather than just the final pacing interval. For instance, in canine ventricles, a pacing-down protocol (with pacing interval decreased in steps of 50 ms to 10 ms) produced a different stationary pattern of alternans than a protocol with constant pacing rate applied to quiescent tissue, for the same asymptotic rate Gizzi et al. (2013). Similar dependence was found in normal rabbit hearts Ziv et al. (2009): while decreasing the pacing interval by 10 ms steps failed to induce discordant alternans, this state was reached by using smaller steps of 5 ms to 2 ms.
The objective of this paper is to establish a connection between transient amplification, memory effects, and the transition to atrial fibrillation in a simple two-dimensional model of cardiac tissue featuring discordant alternans. This paper is structured as follows. Section II describes the model used in this study. The results and their discussion are presented in Sections III and IV, respectively.
Ii Model of paced atrial tissue
Since our focus here is on fundamental mechanisms rather than detailed comparison with experiments, we chose a simple monodomain model of atrial tissue. It has the form of a reaction-diffusion partial differential equation (PDE)
where the field describes the state of cardiac cells (cardiomyocytes), is a diagonal matrix of diffusion coefficients that represents the coupling between adjacent cells, is the ionic model that describes the dynamics of individual cells, and is the (scaled) density of the external pacing current.
Furthermore, we chose the ionic model introduced by Karma Karma (1994), which captures the essential alternans instability responsible for initiating and maintaining complex arrhythmic behaviors. To make it differentiable, the model has been modified Byrne, Marcotte, and Grigoriev (2015); Marcotte and Grigoriev (2015) such that
where , is a scaled transmembrane voltage, and is a gating variable. Unless otherwise stated, the parameters values used were , , , , , , and , where is the restitution parameter. We also used no-flux boundary conditions , where is a unit vector normal to the domain boundary, as a physiologically accurate representation of the boundary of excitable tissue.
At the chosen value of the restitution parameter, in two dimensions, this model reproduces the transition from normal rhythm to alternans to fibrillation observed in experiments as pacing frequency is increased. Yet, with just two variables, as opposed to several tens of variables in the most complex models Fenton and Cherry (2008), the Karma model is simple enough to facilitate the computationally demanding analysis presented in this work.
All numerical simulations were carried out on a square domain of side length cm, which is close to the minimal size required for sustained SWC in Karma model with the present choice of parameters. The PDE (1) was discretized in space using finite differences on a computational grid, which corresponds to the size of one computational cell cm roughly equal to the physical length of one cardiac cell. The values of and on this grid define the state of the system as a point in a 18432-dimensional state space. The discretized equations were solved using 4th order Runge-Kutta method with time step ms chosen to be much smaller than the shortest characteristic time scale of the model, ms.
The pacing was accomplished by applying a current to a patch of grid points located near the upper-left corner of the domain. The patch was displaced by one grid point vertically in order to break the symmetry with respect to reflection about the diagonal of the computational domain. The pacing current was in the form of square pulses of 5 ms duration,
Denoting the time at which the -th pulse is applied, the pacing current can be written as
where for paced grid points and zero otherwise.
iii.1 Wave breakup induced by fast pacing
We investigated stability of the excitation waves produced by pacing with gradually decreasing pacing interval , , starting at . For some tissue models, a decrease in the pacing interval does not produce wave breakup unless (i) the tissue is heterogeneous or (ii) an external current is applied at a location other than the pacing site (ectopic beat) Qu et al. (2000). However, for the Karma model considered here, neither of these conditions is necessary for wave breakup.
Let us define , , to be a sequence of pacing stimuli with a constant interval: . Consider the pacing protocol A: ms, followed by ms, ms, ms, and ms (cf. Fig. 1) with initial condition corresponding to the stable uniform equilibrium . The protocol B is identical to protocol A, except for ms (dashed line in Fig. 1). For protocol B, each pacing stimulus generates a traveling wave with nearly circular wavefront and waveback (cf. Fig. 2(a)) that is absorbed at the bottom and right domain boundaries. The asymptotic state is time-periodic, with period ms, and corresponds to discordant alternans as we will see below. During the first 55 pacing intervals, pacing protocol A produces qualitatively the same dynamics as protocol B. On the 56th pacing interval, however, conduction block causes wave breakup (cf. Fig. 2(b)), creating two spiral waves on the 57th interval (cf. Fig. 2(c)), which undergo additional breakups and eventually transition to SWC (cf. Fig. 2(d)). SWC corresponds to atrial fibrillation; both feature re-entrant waves and are sustained irrespective of the pacing.
So how can we understand the difference in the outcomes of the two protocols? Although this difference is a clear indication that the memory effect plays an important role in our system, this does not explain much. The memory is actually quite long: indeed, 15(!) pacing stimuli preceding conduction block in protocol A were delivered at the same asymptotic rate. Given that the asymptotic alternans state is stable, it would have been natural to expect that any deviation present at would all but disappear after such a long period. In fact, as we show below, infinitesimal initial perturbations about this stable alternans state can grow, albeit transiently, by several orders of magnitude over long time intervals. Small, but finite, initial perturbations can grow enough to cause a transition to a completely different (here chaotic) attractor, as it happens for protocol . This is a signature of multistability often associated with subcritical instabilities; similar phenomenology is found during bypass transition to turbulence in fluid flows Orszag and Patera (1980); Chaté and Manneville (1987).
iii.2 Periodic solutions
Let us start the analysis of the problem with a description of the various solutions to (1). To explicitly express the dependence of a solution of (1) on the initial condition , we define the time evolution operator
For constant pacing, , various time-periodic solutions can be defined as fixed points or -cycles of the Poincaré map , which describes the discrete-time dynamics on the Poincaré section . States with different periods corresponds to solutions of the nonlinear equation
with different . In particular, the 1:1 response corresponds to the fixed point of the Poincaré map . At faster pacing rates, two different solutions of (6) with appear. The 2-cycle , where and , corresponds to the alternans solution (2:2 response). A different 2-cycle , where and , corresponds to the conduction block solution (2:1 response). These three periodic solutions are compared in Fig. 3.
Stability of time-periodic solutions can be described in terms of the dynamics of deviations from the reference solution . The dynamics of finite disturbances is described by
For infinitesimal disturbances, (7) reduces to a linear PDE
where and is the time-dependent Jacobian of the nonlinear term describing cellular dynamics evaluated on the periodic orbit.
The linear time-evolution operator that advances the solution of (8) from to in the tangent space,
can be found by integrating (8), which yields a time-ordered exponential
Linear stability of is determined by the eigenvalues (Floquet multipliers) of the operator ,
where are the corresponding eigenfunctions (Floquet modes). In particular, if for all , any infinitesimal disturbance will eventually decay to zero, so that for . In the following discussion, we will index the eigenvalues in the order of decreasing absolute value, .
Asymptotic stability, however, does not imply that initial disturbances will decay monotonically. If the evolution operator is nonnormal (i.e., ), some disturbances may grow transiently, before asymptotic decay takes over Trefethen et al. (1993). As we have shown previously Garzón, Grigoriev, and Fenton (2011, 2014) in the context of one-dimensional paced cardiac tissue (Purkinje fibers), such transient growth may arise, for example, due to the action of closed-loop feedback control, which makes the 1:1 solution linearly stable. In the present case, the Jacobian is non-self-adjoint, so both the differential operator and the evolution operator are nonnormal, hence, generalized linear stability theory Ioannou and Farrell (1996a, b) is required to describe transient response to perturbations associated with changes in the pacing rate.
The adjoint of the evolution operator
plays an important role in the generalized stability theory. It defines the evolution of disturbances in the adjoint of the tangent space,
backward in time. In practice it is more straightforward to compute the action of on a vector by solving the linear PDE
with “initial” condition defined at the final time.
Given that time-periodic states could be stable or unstable, we computed them as solutions of (6) using a Newton-Krylov solver Garzón, Grigoriev, and Fenton (2011); Marcotte and Grigoriev (2015). The spectrum of the evolution operator was computed using the Arnoldi method implemented by the Matlab function eigs Lehoucq and Sorensen (1996) with accuracy yielding eight significant digits. In both instances matrix-free evaluation of the the matrix-vector product was accomplished via simultaneous numerical time-integration of (1) and (8). The matrix-free calculation of the product involves backward integration and cannot be accomplished by integrating (1) and (14) simultaneously. Instead a precomputed and interpolated solution was used to compute the solution to (14) using the procedure described in Ref. Marcotte and Grigoriev, 2016. The time integrators were implemented as OpenCL kernels Marcotte and Grigoriev (2012) administered by host codes written in C and encapsulated in Matlab mex-files. This allowed substantially speeding up the calculations by executing them on general purpose Graphics Processing Units (GPUs), at the same time retaining convenient interface within Matlab scripts.
iii.3 Stability spectra
The 1:1 solution exists at all considered in this study. It is stable for ms and becomes unstable for . The leading eigenvalue crosses the unit circle along the negative real line, which corresponds to a period doubling bifurcation, at which point the alternans solution with period is created. It corresponds to discordant alternans with four stationary nodal lines, as Fig. 4(a) illustrates. The alternans is created discordant, since the domain is sufficiently large, according to theoretical analysis Echebarria and Karma (2002a, 2007). In experiments Pastore et al. (1999); Walker and Rosenbaum (2003) and simulations Qu et al. (2000); Watanabe et al. (2001); Gizzi et al. (2013) alternans is typically created concordant and becomes discordant as the pacing interval is decreased. However, this subtle distinction has no bearing on the development of conduction block that plays the crucial role in transition to fibrillation.
To further characterize this bifurcation, we computed the distance between the 1:1 and 2:2 solutions in the state space, defined using the Euclidean 2-norm
where the scalar product is given by
and the integral is over the entire computational domain . The bifurcation is supercritical (cf. Fig. 4(b)), also in agreement with the theoretical predictionsEchebarria and Karma (2002a), and the 2:2 solution exists for all (at least down to ms).
The stability balloon for the 2:2 solution in the plane is shown in Fig. 5. This solution undergoes a subcritical Hopf bifurcation when the neutral stability curve (shown as a solid line) is crossed. The minimum of the neutral stability curve is at and ms. The value considered here is below , so the 2:2 solution turns out to be linearly stable for all . The corresponding Floquet spectra for different values of are shown in Fig. 6. The change in the leading eigenvalue associated with a decrease in is shown in Fig. 6(a). Although this eigenvalue never leaves the unit circle, it approaches this circle for ms.
Above the neutral stability curve infinitesimally small perturbations of the 2:2 state cause a direct transition to SWC. This chaotic state exists over a region in the plane that is bigger than the gray-shaded area in Fig. 5, so our system indeed features both a subcritical instability and multistability – the hallmarks of a bypass transition. Of particular relevance for this study, persistent SWC exists for and ms. In addition, there is one more stable state that exists for the same parameter values, the one that describes a 2:1 response. For , it is created in a saddle-node bifurcation at ms and remains stable as decreases to at least 50 ms (that is, over the entire range of in Fig. 5).
Given this information, how can we understand the difference in the outcomes of the pacing protocols A and B? Both protocols can be viewed as being composed of two stages: the first stage (first 40/50 pacing intervals for protocol A/B) determines the initial condition for the second stage (pacing with constant ms that starts, respectively at or ). The perturbation (cf. Fig. 7(c-d)) in protocol B eventually decays, in agreement with the prediction of linear stability analysis, while the perturbation (cf. Fig. 7(a-b)) in protocol A grows, producing conduction block and wave breakup after 15 pacing intervals and, eventually, transition to persistent SWC. (Note that in the 2:1 state conduction block is global, annihilating an entire excitation wave, while in protocol A conduction block is local and leads to a breakup, but not a total annihilation, of an excitation wave.)
Since three different stable solutions (2:2, 2:1, and SWC) coexist at ms, the asymptotic regime is determined by the initial condition. In particular, there is a neighborhood of the 2-cycle , such that for all initial conditions , the orbit approaches this 2-cycle as . is known as the basin of attraction of the 2-cycle and has a finite size. Similarly, the 2-cycle has a basin of attraction . In particular, the initial condition in protocol B lies inside , so the asymptotic state corresponds to the 2:2 solution. On the other hand, the initial condition in protocol A lies outside of both and , so the asymptotic state corresponds to SWC.
Since the evolution operator describing the linearized dynamics is nonnormal, the basin of attraction becomes quite thin in some directions. Fig. 8 shows a two-dimensional cartoon of this; the actual dimensionality of is the same as that of the state space of the discretized PDE, 18432 in the present case. We are primarily interested in the direction in the state space (which corresponds to perturbation shape in the physical space) for which the distance from to the boundary of the basin of attraction is the smallest (we will refer to it as the optimal direction), since the dynamics are the most sensitive to perturbations corresponding to this direction. The point on that is the closest to the 2-cycle defines the critical optimal perturbation . Although the shape of , and hence the perturbation , is defined by the nonlinear evolution operator (5), as we illustrate below, both the direction of and the strong transient growth associated with perturbations in this direction can be understood reasonably well within the linear approximation.
iii.4 Generalized linear stability analysis
Since is self-adjoint, it has real eigenvalues and orthogonal eigenvectors ,
with eigenvalues sorted from largest to smallest. Expanding the initial condition in the eigenvector basis yields
where , such that
The maximum of over all initial disturbances with a fixed norm is reached when all , except for , i.e., .
Let and define . The action of the evolution operator can now be expressed in a simple manner. Multiplying both sides of (19) by and recognizing that is arbitrary, we obtain
which is the singular value decomposition (SVD) of with , , and being, respectively, the singular values, left singular vectors, and right singular vectors of . The leading singular vectors and singular value have a useful interpretation: as we have already found, the right singular vector determines the shape of the optimal initial perturbation that is amplified the most (in the sense of the 2-norm (17)) at time . The left singular vector determines the shape this optimal perturbation acquires at time . Finally, the singular value determines the magnitude of transient amplification which corresponds to this particular initial condition.
For small initial disturbances, the time at which the maximal transient amplification is achieved is determined by the balance between transient growth and asymptotic decay. The weaker the asymptotic decay, the larger the amplification that can be achieved. Computation of is numerically rather expensive, since SVD requires alternate direction integration of the linear PDEs (8) and (14) on the time interval and there is no simple relation that allows one to compute the singular values at time based on their values at time . We therefore performed this calculation over many pacing periods on a coarse grid with .
The leading singular value is shown in Fig. 9. We see that transient growth is a relatively slow process. The maximal transient amplification takes place only after almost 30 pacing stimuli. At this level of resolution, appears to be a discontinuous function of time. In fact, this is not the case, as the time-resolved calculation over the interval illustrates (cf. Fig. 10). We find that has a rather fine structure on a time scale smaller than one pacing period. The maxima of transient amplification have a pattern that is -periodic; they are located at , where is an integer. The peaks’ values, on the other hand, are modulated with a period determined by the leading Floquet multiplier. The largest transient amplification () is achieved at and the second highest value is achieved at . Recall that for protocol A it takes a time interval of roughly until the dynamics transitions from the neighborhood of the 2:2 state to SWC. These times are all of the same order of magnitude, suggesting that memory effect extends over times scales of order ms, much longer than the naive estimate ms would suggest.
The right singular vectors (cf. Fig. 11(a-b)) associated with the time instances when is near the peak value are indistinguishable, implying that their spatial structure, which determines the optimal initial perturbation is a robust feature of the dynamics. Indeed, the time-dependence of for can be reproduced with extremely high accuracy by the norm of the disturbance which corresponds to the optimal one, e.g., , where is the amplitude of the initial disturbance. As Fig. 10 illustrates, the norm is indistinguishable from once both have been normalized by their maximal values. Effectively, this means that the initial disturbance is in fact optimal for all times . This optimal disturbance is strongly localized around the pacing site, which means that the dynamics are very sensitive to the perturbation in the timing of the pacing stimuli.
The corresponding left singular vector is shown in Fig. 11(c-d)). It is worth noting that, unlike the right singular vector which is strongly localized, the left singular vector is not. This spatial delocalization is to a large extent responsible for the apparent large transient amplification with respect to the 2-norm (15), which involves integration over the entire spatial domain. Furthermore, the right singular vector is dominated by the component, so that the dynamics are more sensitive to the perturbation of the variable than the variable. For the left singular vector the opposite is true, hence the perturbation affects the variable to a much larger degree than the variable. A similar asymmetry was found for the left and right marginal eigenvectors (response functions and Goldstone modes) of spiral wave solutions, which similarly describe the cause-and-effect relationship for small perturbations Marcotte and Grigoriev (2016).
As expected, generalized linear stability analysis correctly predicts the dynamics produced by optimal initial disturbances of sufficiently small amplitude . Evolving the optimal initial perturbation with magnitude using the nonlinear equation (7), we find that grows transiently, achieving the predicted amplification at time (cf. Fig. 12(a)). For the asymptotic exponential decay with the predicted rate sets in, where is the leading Floquet multiplier. A stronger perturbation with leads to the same result, with the entire curve simply shifted up, indicating that nonlinear effects are still negligible.
These results are qualitatively similar to transient amplification found for traveling wave solutions describing thin fluid films spreading on a solid substrate Bertozzi and Brenner (1997). In that case as well, transient growth of small localized perturbations was traced to the nonnormality of the evolution operator Grigoriev (2005). Strong transient amplification of disturbances also characterizes traveling excitation wave in models of paced Purkinje fibers with feedback, although in that case the base state corresponds to the 1:1 response Garzón, Grigoriev, and Fenton (2011, 2014). In fact, increasing transient amplification was identified as the reason for breakdown of feedback control for long fibers controlled at the pacing site.
iii.5 Finite amplitude disturbances
Since the 2:2 state is linearly stable, transition to SWC is subcritical and requires a finite amplitude perturbation. The critical magnitude of the perturbation away from that causes a transition to a different solution (, , , or SWC) is given by the distance from to the nearest boundary of the corresponding basin of attraction (cf. Fig. 8). For linearly optimal initial perturbations, , we find that the asymptotic state is for , as Fig. 12(b) shows. For the asymptotic state is , which is also a 2:2 state, but with a different phase. For the asymptotic state is again , but for the asymptotic state is SWC, which persists for at least . Increasing further we find either persistent SWC or transient SWC that eventually gives way to one of the 2:1 states ( or ).
The extreme sensitivity of the outcome to the choice of initial conditions is a characteristic feature of chaotic attractors that have fractal basin boundaries McDonald et al. (1985). Since and are parts of the same 2-cycle, both and share a part of their boundary with the basin of attraction of the SWC state. Like a quilt, basin boundaries are composed of a patchwork of pieces, each one corresponding to the stable manifold of an edge state – a saddle solution embedded into the boundary, which has a one-dimensional unstable manifold. The relation between optimal perturbations, transient growth, edge states, and basin boundaries has originally been investigated in the context of simple nonlinear PDEs Handel and Grigoriev (2006). More recently, edge states have been studied extensively as gateways mediating subcritical transition to turbulence in shear fluid flows Schneider, Eckhardt, and Yorke (2007); Gibson, Halcrow, and Cvitanović (2008).
The 1:1 state is an example of such an edge state. It has one unstable direction and its stable manifold forms the boundary between and . As Fig. 12(b) illustrates, the solutions corresponding to perturbations with and approach and closely follow the same periodic solution for , before eventually separating and approaching, respectively, and . This periodic solution is indeed , as shown in Fig. 13. As expected, the two trajectories approach along the least stable direction and separate along the unstable direction, with the rates predicted by the linearization.
The basin boundary of the chaotic attractor lies further away from than the basin boundary of the 2:2 state . An upper bound for the distance to is given by the norm of the critical linearly optimal initial disturbance. It is worth noting that the norm of this disturbance is more than two orders of magnitude less than the norm of the state itself. The ratio of the norms is not as large as the maximal transient amplification , but it is of the same order of magnitude, which suggests that linear theory cold be used to estimate the critical disturbance magnitude. However, a fully nonlinear calculation is required to obtain a quantitatively accurate prediction for both the magnitude and the shape of the smallest initial disturbance required to trigger a transition from the 2:2 state to SWC.
Figure 14 shows the evolution of the initial condition , which corresponds to the linearly optimal initial disturbance with a magnitude that places it just outside of . Even though the 2:2 state is linearly stable and the perturbation is quite small, transient amplification of the initial disturbance leads to pronounced alternans which causes conduction block in the lower right corner of the domain at . Recall that, according to generalized linear stability theory, the shape of the perturbation at the moment when it reaches the largest amplification is given by the left singular vector of : . The location of conduction block is indeed found to be consistent with the shape of the left singular vector , which 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 (cf. Fig. 15) is noticeably different from (cf. Fig. 11(c,d)).
While the first instance of conduction block generates a pair of spiral waves at , this does not cause an immediate transition to SWC. The spiral waves have too little room to sustain themselves: they collide with the next wavefront emitted by the pacing site and are annihilated at . Another instance of conduction block occurs at , again generating two spiral waves at . Now conduction block occurs sooner in the cycle and further away from the lower right corner of the domain, and the spiral waves survive the collision with the next wavefront emitted by the pacing site. (Note, however, that in both instances conduction block occurs near the nodal line that is the furtherest from the pacing site.) Since the period of the spiral waves is shorter than the pacing period, they gradually invade the domain, generating further breakups and initiating sustained SWC.
As we discussed previously, from on, protocol A corresponds to constant pacing with ms. The initial condition for the last stage of protocol A corresponds to a fairly significant perturbation away from the 2:2 state . The norm of this perturbation is two orders of magnitude larger than (and is close to the norm of itself). On the other hand, the spatial structure of this perturbation (cf. Fig. 7(a-b)) is quite different from that of the optimal perturbation (cf. Fig. 11(a-b)) and therefore it does not experience transient amplification, as Fig. 16 illustrates. This difference can be quantified using the scalar product , where the hat denotes normalization (e.g., ), which shows that has essentially no component along . At the same time , so that has a significant component in the direction of (cf. Fig. 8).
Despite the dramatic difference in the amplitude and spatial structure of the perturbations and , 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 possible explanation is that, before approaching the chaotic attractor, both trajectories closely follow the boundary between and for an extended period of time (around for both and ) and finally they start to separate from it along the unstable manifold of the same edge state that lies on this boundary.
Although itself unstable, the basin boundary can be followed using the bisection algorithm Schneider, Eckhardt, and Yorke (2007) that has been successfully used in the context of fluid turbulence. For the system considered here we found that the critical optimal perturbation approaches a chaotic edge state that lies on the boundary of . It may in principle be possible to find nonchaotic edge states embedded into as well. For instance, the unstable quasiperiodic solution born an the subcritical Hopf bifurcation of the 2:2 state lies on this boundary and, if it is has one unstable eigenvalue, it would be an example of an edge state. Even simpler (e.g., time-periodic) edge states may be found by applying the bisection algorithm to initial conditions other than , however this is far outside the scope of this study.
The perturbation that corresponds to protocol B is also much stronger than the critical optimal one, but since its spatial structure (cf. Fig. 7(c-d)) is also quite different from , it does not experience transient amplification either, as Fig. 16 illustrates. Since and , has essentially no component along , but a significant component along . The norm of the perturbation, , is small enough for the initial condition to fall within the basin of attraction of the 2:2 state (cf. Fig. 8). Hence the perturbation eventually decays, and the dynamics return to concordant alternans instead of transitioning to SWC.
The results presented in this paper help improve our understanding of the transition to fibrillation in paced atrial tissue and the role of alternans. Specifically, these results confirm that the first step in the transition caused by an increase in the pacing rate is the development of discordant alternans Pastore et al. (1999); Weiss et al. (2006); Gizzi et al. (2013). Theoretical analysis of alternans in one dimension predicts Echebarria and Karma (2002a); Fox, Gilmour Jr, and Bodenschatz (2002) that conduction block should occur at the location that is the furthest from the pacing site. On the other hand, experimental and numerical evidence suggests that reentry is promoted by steep repolarization gradients Konta et al. (1990); Tachibana et al. (1998); Weiss et al. (2006) and therefore conduction block should happen near a nodal line. Indeed, we find that conduction block that initiates wave breakup and reentry happens to lie close to the nodal line of the stable discordant alternans that is the furtherest from the pacing site. It should be pointed out, however, that this correspondence is not precise, since nodal lines can (and do) move during the transient evolution from stable discordant alternans to conduction block.
However, in a marked departure from the prevailing paradigm, the transition from discordant alternans to conduction block (followed by wave breakup, reentry, and SWC) is not associated with a bifurcation leading to either an instability or disappearance of the alternans state. Instead the transition is triggered by strong transient amplification of disturbances away from stable alternans. At the pacing interval of 80 ms, certain infinitesimal disturbances were found to be amplified by almost three orders of magnitude. Transient amplification increases as the pacing interval is decreased, reaching a maximum of four orders of magnitude for a pacing interval of 66 ms. This is very similar to the trend found for linearly stable fluid flows Meseguer and Trefethen (2003), where transient amplification increases with the Reynolds number, which plays the role analogous to the pacing rate.
Furthermore, we found that the spatial profile of the optimal perturbations that experience the strongest transient amplification corresponds to a tiny change in the pacing interval. Hence even small changes in the pacing interval will generate a strong perturbation away from the 2:2 response. For a protocol with a decreasing pacing interval, the strength of a perturbation is controlled by the rate at which the pacing interval changes. Since transient amplification increases with decreasing pacing interval, transition from alternans to SWC will happen sooner (for a longer pacing interval) if the pacing interval is decreased faster.
Although generalized linear stability analysis is capable of describing transient dynamics qualitatively, it cannot produce a quantitatively accurate description of the transition, which involves intrinsically nonlinear effects, such as conduction block. Using numerical simulations, we were able to predict that, for a pacing interval of 80 ms, transition to SWC is triggered by a linearly optimal perturbation the norm of which is two orders of magnitude smaller than the norm of the alternans state itself. Even smaller, nonlinearly optimal perturbations about the alternans state could trigger a transition to SWC. Computing their magnitude and spatial structure, however, would require a fully nonlinear analysis, following, for example, an adjoint-based minimization approach used to characterize transition to turbulence in shear fluid flows Pringle and Kerswell (2010).
Our results also have interesting implications for the problem of low-energy defibrillation. Although some of the modern approaches Winkle et al. (1989); Fenton et al. (2009); Luther et al. (2011) have demonstrated that the energy required for defibrillation could be reduced substantially compared with the standard defibrillation protocols, further dramatic reduction may be possible by exploiting the geometric structure of the state space. Understanding the dynamics in the region of state space surrounding an edge state that lies on the basin boundary between SWC and time-periodic solutions can help design protocols that can terminate fibrillation with, technically, infinitely small perturbations. In practice, the perturbation strength will depend on how close the chaotic dynamics approach the basin boundary, but still there is a potential for orders-of-magnitude improvement.
The identification of edge states should also help predict the opposite process – when a seemingly tiny perturbation to a stable time-periodic rhythm leads to a transition into SWC. More immediately, our results also support an approach designed to prevent fibrillation, rather than terminate fibrillation after it started. In that approach feedback control is used to stabilize the normal rhythm, thereby preventing the development of alternans in the first place Echebarria and Karma (2002b); Dubljevic, Lin, and Christofides (2008); Garzón, Grigoriev, and Fenton (2011, 2014). This approach should work regardless of whether the alternans state is stable or not (e.g., for in the model considered here).
The analysis presented here was based on a greatly simplified model of cardiac tissue, so while its results may have a profound impact on how we view the dynamical mechanisms that describe the onset of fibrillation, they need to be verified using more detailed electrophysiological models of cardiac tissue in two and three dimensions. In particular, it is important to check whether our main results can be reproduced in bi-domain models that are required to correctly describe intra- and extra-cellular potentials. Finally, it would be interesting to see how the difference between atria and ventricles (both in terms of the electrophysiology and in terms of dimensionality) affect the dynamical description of transition from discordant alternans to fibrillation.
Acknowledgements.This material is based upon work supported by the National Science Foundation under Grant No. CMMI-1028133. The Tesla K20 GPUs used for this research were donated by the “NVIDIA Corporation” through the academic hardware donation program.
- Mozaffarian et al. (2016) D. Mozaffarian, E. J. Benjamin, A. S. Go, D. K. Arnett, et al., “Heart disease and stroke statistics-2016 update a report from the American Heart Association,” Circulation 133, e38–e360 (2016).
- Chugh et al. (2014) S. S. Chugh, R. Havmoeller, K. Narayanan, D. Singh, et al., “Worldwide epidemiology of atrial fibrillation: A global burden of disease 2010 study,” Circulation 129, 837–847 (2014).
- Nolasco and Dahlen (1968) J. B. Nolasco and R. W. Dahlen, ‘‘A graphic method for the study of alternation in cardiac action potentials,” J. Appl. Physiol. 25, 191–196 (1968).
- Rosenbaum et al. (1994) D. S. Rosenbaum, L. E. Jackson, J. M. Smith, H. Garam, J. N. Ruskin, and R. J. Cohen, “Electrical alternans and vulnerability to ventricular arrythmias,” N. Engl. J. Med. 330, 235 (1994).
- Bloomfield et al. (2006) 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 alternans and the risk of death or sustained ventricular arrhythmias in patients with left ventricular dysfunction,” Journal of the American College of Cardiology 47, 456–463 (2006).
- Hiromoto et al. (2005) K. Hiromoto, H. Shimizu, Y. Furukawa, T. Kanemori, T. Mine, T. Masuyama, and M. Ohyanagi, “Discordant repolarization alternans-induced atrial fibrillation is suppressed by verapamil,” Circ J 69, 1368–1373 (2005).
- Narayan et al. (2011) S. M. Narayan, M. R. Franz, P. Clopton, E. J. Pruvot, and D. E. Krummen, “Repolarization alternans reveals vulnerability to human atrial fibrillation,” Circulation 123, 2922–2930 (2011).
- Narayan and Zaman (2016) S. Narayan and J. Zaman, “Mechanistically based mapping of human cardiac fibrillation,” J Physiol 594, 2399–2415 (2016).
- Majumder et al. (2016) R. Majumder, M. C. Engels, A. A. F. de Vries, A. V. Panfilov, and D. A. Pijnappels, “Islands of spatially discordant APD alternans underlie arrhythmogenesis by promoting electrotonic dyssynchrony in models of fibrotic rat ventricular myocardium,” Scientific Reports 6, 24334 (2016).
- Chang and Trayanova (2016) K. C. Chang and N. A. Trayanova, “Mechanisms of arrhythmogenesis related to calcium-driven alternans in a model of human atrial fibrillation,” Scientific Reports 6, 36395 (2016).
- Pastore et al. (1999) J. M. Pastore, S. D. Girouard, K. R. Laurita, F. G. Akar, and D. S. Rosenbaum, “Mechanism linking T-wave alternans to the genesis of cardiac fibrillation,” Circulation 99, 1385–1394 (1999).
- Weiss et al. (2006) J. N. Weiss, A. Karma, Y. Shiferaw, P. S. Chen, A. Garfinkel, and Z. Qu, “From pulsus to pulseless: The saga of cardiac alternans,” Circulation Research 98, 1244–1253 (2006).
- Echebarria and Karma (2002a) B. Echebarria and A. Karma, “Instability and spatiotemporal dynamics of alternans in paced cardiac tissue,” Phys. Rev. Lett. 88, 208101 (2002a).
- Echebarria and Karma (2007) B. Echebarria and A. Karma, “Amplitude equation approach to spatiotemporal dynamics of cardiac alternans,” Phys. Rev. E 76, 051911 (2007).
- Skardal, Karma, and Restrepo (2014) P. S. Skardal, A. Karma, and J. G. Restrepo, “Spatiotemporal dynamics of calcium-driven cardiac alternans,” Physical Review E 89, 052707 (2014).
- Fenton et al. (2002) F. H. Fenton, E. M. Cherry, H. M. Hastings, and S. J. Evans, “Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity,” Chaos 12, 852–892 (2002).
- Karma (1994) A. Karma, “Electrical alternans and spiral wave breakup in cardiac tissue,” Chaos 4, 461–472 (1994).
- Fox, Gilmour Jr, and Bodenschatz (2002) J. J. Fox, R. F. Gilmour Jr, and E. Bodenschatz, “Conduction block in one-dimensional heart fibers,” Physical Review Letters 89, 198101 (2002).
- Alonso, Bär, and Echebarria (2016) S. Alonso, M. Bär, and B. Echebarria, “Nonlinear physics of electrical wave propagation in the heart: a review,” Reports on Progress in Physics 79, 096601 (2016).
- Surovyatkina et al. (2010) E. Surovyatkina, D. Noble, D. Gavaghan, and A. Sher, “Multistability property in cardiac ionic models of mammalian and human ventricular cells,” Progress in Biophysics and Molecular Biology 103, 131–141 (2010).
- Skardal and Restrepo (2014) P. S. Skardal and J. G. Restrepo, “Coexisting chaotic and multi-periodic dynamics in a model of cardiac alternans,” Chaos 24, 043126 (2014).
- Comtois and Vinet (2005) P. Comtois and A. Vinet, “Multistability of reentrant rhythms in an ionic model of a two-dimensional annulus of cardiac tissue,” Physical Review E 72, 051927 (2005).
- Oliver, Henriquez, and Krassowska (2005) R. A. Oliver, C. S. Henriquez, and W. Krassowska, “Bistability and correlation with arrhythmogenesis in a model of the right atrium,” Annals of Biomedical Engineering 33, 577–589 (2005).
- Henningson, Lundbladh, and Johansson (1993) D. S. Henningson, A. Lundbladh, and A. V. Johansson, “A mechanism for bypass transition from localized disturbances in wall-bounded shear flows,” Journal of Fluid Mechanics 250, 169–207 (1993).
- Trefethen et al. (1993) L. N. Trefethen, A. E. Trefethen, S. Reddy, and T. A. Driscoll, “Hydrodynamic stability without eigenvalues,” Science 261, 578–584 (1993).
- Chialvo, Michaels, and Jalife (1990) D. R. Chialvo, D. C. Michaels, and J. Jalife, “Supernormal excitability as a mechanism of chaotic dynamics of activation in cardiac purkinje fibers.” Circulation Research 66, 525–545 (1990).
- Rosenbaum et al. (1982) M. B. Rosenbaum, H. H. Blanco, M. V. Elizari, J. O. Lázzari, and J. M. Davidenko, “Electrotonic modulation of the t wave and cardiac memory,” Am. J. Cardiol. 50, 213–222 (1982).
- Hund and Rudy (2000) T. J. Hund and Y. Rudy, “Determinants of excitability in cardiac myocytes: mechanistic investigation of memory effect,” Biophys. J. 79, 3095–3104 (2000).
- Gizzi et al. (2013) A. Gizzi, E. M. Cherry, R. F. Gilmour, S. Luther, S. Filippi, and F. H. Fenton, “Effects of pacing site and stimulation history on alternans dynamics and the development of complex spatiotemporal patterns in cardiac tissue,” Frontiers in Physiology 4, 71 (2013).
- Ziv et al. (2009) 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 complex behaviour of spatially discordant alternans in a transgenic rabbit model of type 2 long QT syndrome.” The Journal of physiology 587, 4661–4680 (2009).
- Byrne, Marcotte, and Grigoriev (2015) G. Byrne, C. D. Marcotte, and R. O. Grigoriev, “Exact coherent structures and chaotic dynamics in a model of cardiac tissue,” Chaos 25, 033108 (2015).
- Marcotte and Grigoriev (2015) C. D. Marcotte and R. O. Grigoriev, “Unstable spiral waves and local Euclidean symmetry in a model of cardiac tissue,” Chaos 25, 063116 (2015).
- Fenton and Cherry (2008) F. H. Fenton and E. M. Cherry, “Models of cardiac cell,” Scholarpedia 3, 1868 (2008).
- Qu et al. (2000) Z. Qu, A. Garfinkel, P. S. Chen, and J. N. Weiss, “Mechanisms of discordant alternans and induction of reentry in simulated cardiac tissue,” Circulation 102, 1664–1670 (2000).
- Orszag and Patera (1980) S. A. Orszag and A. T. Patera, “Subcritical transition to turbulence in plane channel flows,” Phy. Rev. Lett. 45, 989 (1980).
- Chaté and Manneville (1987) H. Chaté and P. Manneville, “Transition to turbulence via spatio-temporal intermittency,” Phys. Rev. Lett. 58, 112 (1987).
- Garzón, Grigoriev, and Fenton (2011) A. Garzón, R. O. Grigoriev, and F. H. Fenton, “Model-based control of cardiac alternans in Purkinje fibers,” Phys. Rev. E 84, 041927 (2011).
- Garzón, Grigoriev, and Fenton (2014) A. Garzón, R. O. Grigoriev, and F. H. Fenton, “Continuous-time control of alternans in long purkinje fibers,” Chaos 24, 033124 (2014).
- Ioannou and Farrell (1996a) P. J. Ioannou and B. F. Farrell, “Generalized stability theory part 1: Autonomous operators,” J. Atmos. Sci. 53, 2025–2040 (1996a).
- Ioannou and Farrell (1996b) P. J. Ioannou and B. F. Farrell, “Generalized stability theory part 2: Nonautonomous operators,” J. Atmos. Sci. 53, 2041–2053 (1996b).
- Lehoucq and Sorensen (1996) R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly re-started Arnoldi iteration,” SIAM J. Matr. Anal. Appl. 17, 789 (1996).
- Marcotte and Grigoriev (2016) C. D. Marcotte and R. O. Grigoriev, “Adjoint eigenfunctions of temporally recurrent single-spiral solutions in a simple model of atrial fibrillation,” Chaos , 093107 (2016).
- Marcotte and Grigoriev (2012) C. Marcotte and R. O. Grigoriev, “Implementation of PDE models of cardiac dynamics on GPUs using OpenCL,” (2012), arXiv:1309.1720.
- Walker and Rosenbaum (2003) M. L. Walker and D. S. Rosenbaum, “Repolarization alternans: implications for the mechanism and prevention of sudden cardiac death,” Cardiovasc. Res. 57, 599 (2003).
- Watanabe et al. (2001) M. A. Watanabe, F. H. Fenton, S. J. Evans, H. M. Hastings, and A. Karma, “Mechanisms for discordant alternans,” J. Cardiovasc. Electrophysiol. 12, 196–206 (2001).
- Bertozzi and Brenner (1997) A. L. Bertozzi and M. P. Brenner, “Linear stability and transient growth in driven contact lines,” Physics of Fluids 9, 530–539 (1997).
- Grigoriev (2005) R. O. Grigoriev, “Transient growth in driven contact lines,” Physica D 205, 105–116 (2005).
- McDonald et al. (1985) S. W. McDonald, C. Grebogi, E. Ott, and J. A. Yorke, “Fractal basin boundaries,” Physica D 17, 125–153 (1985).
- Handel and Grigoriev (2006) A. Handel and R. O. Grigoriev, “Transient dynamics and nonlinear stability of spatially extended systems,” Phys. Rev. E 74, 036302 (2006).
- Schneider, Eckhardt, and Yorke (2007) T. M. Schneider, B. Eckhardt, and J. A. Yorke, “Turbulence transition and the edge of chaos in pipe flow,” Phys. Rev. Lett. 99, 034502 (2007).
- Gibson, Halcrow, and Cvitanović (2008) J. F. Gibson, J. Halcrow, and P. Cvitanović, “Visualizing the geometry of state-space in plane Couette flow,” J. Fluid Mech. 611, 107–130 (2008).
- Konta et al. (1990) T. Konta, K. Ikeda, M. Yamaki, K. Nakamura, K. Honma, I. Kubota, and S. Yasui, “Significance of discordant st alternans in ventricular fibrillation.” Circulation 82, 2185–2189 (1990).
- Tachibana et al. (1998) H. Tachibana, I. Kubota, M. Yamaki, T. Watanabe, and H. Tomoike, “Discordant s-t alternans contributes to formation of reentry: a possible mechanism of reperfusion arrhythmia,” Am. J. Physiol. 275, H116 (1998).
- Meseguer and Trefethen (2003) A. Meseguer and L. N. Trefethen, “Linearized pipe flow to reynolds number ,” Journal of Computational Physics 186, 178–197 (2003).
- Pringle and Kerswell (2010) C. C. T. Pringle and R. R. Kerswell, “Using nonlinear transient growth to construct the minimal seed for shear flow turbulence,” Phys. Rev. Lett. 105, 154502 (2010).
- Winkle et al. (1989) R. A. Winkle, R. H. Mead, M. A. Ruder, V. Gaudiani, W. S. Buch, B. Pless, M. Sweeney, and P. Schmidt, “Improved low energy defibrillation efficacy in man with the use of a biphasic truncated exponential waveform,” American Heart Journal 117, 122–127 (1989).
- Fenton et al. (2009) F. H. Fenton, S. Luther, E. M. Cherry, N. F. Otani, V. Krinksy, A. Pumir, E. Bodenschatz, and R. F. G. Jr., “Termination of atrial fibrillation using pulsed low-energy far-field stimulation,” Circulation 120, 467–476 (2009).
- Luther et al. (2011) 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 electrical turbulence in the heart,” Nature 475, 235–9 (2011).
- Echebarria and Karma (2002b) B. Echebarria and A. Karma, “Spatiotemporal control of cardiac alternans,” Chaos 12, 923–930 (2002b).
- Dubljevic, Lin, and Christofides (2008) S. Dubljevic, S.-F. Lin, and P. D. Christofides, “Studies of feedback control of cardiac alternans,” Comp. Chem. Eng. 32, 2086 (2008).