Memory effects, transient growth, and wave breakup in a model of paced atrium

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,¹ or the atria, where it creates significant health complications.² 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,³ 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^4,5 and in the atria.^6–8 Strong correlation between alternans and fibrillation was also found in the simulations employing high-fidelity models of two-dimensional ventricular⁹ and atrial¹⁰ tissue.

In tissue, alternans manifests as spatial modulation in the width of the excitation waves,¹¹ 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.¹² Theoretical analysis of propagating excitation waves in one dimension based on amplitude equations^13,14 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 for calcium-driven alternans,¹⁵ which qualitatively behave like their voltage-driven cousins. Theoretical studies of alternans in higher dimensions¹⁴ are too limited in scope to make any specific predictions regarding the transition from alternans to fibrillation, but based on experiments and numerical simulations¹⁶ we know that conduction block that leads to wave breakup plays a key role.

As the pacing rate or the tissue size is increased, conduction block is predicted to arise when alternans becomes unstable, leading to unbounded growth of the amplitude of modulation.^13,17,18 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).¹⁹ 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,²⁰ one,²¹ and two^22,23 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²⁴ and involves strong transient amplification of small disturbances.²⁵ 

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.^26–28 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 a pacing interval decreased in steps of either 50 ms or 10 ms) produced a different stationary pattern of alternans than a protocol with a constant pacing rate applied 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 failed to induce discordant alternans, this state was reached by using smaller steps of either 5 ms or 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 Secs. III and IV, respectively.

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 $u=[u,v](r,t)$ describes the state of cardiac cells (cardiomyocytes), D is a diagonal matrix of diffusion coefficients that represents the coupling between adjacent cells, $f(u)$ is the ionic model that describes the dynamics of individual cells, and I_p is the (scaled) density of the external pacing current.

Furthermore, we chose the ionic model introduced by Karma,¹⁷ which captures the essential alternans instability responsible for initiating and maintaining complex arrhythmic behaviors. To make it differentiable, the model has been modified^31,32 such that

where $Θα(u)=[1+tanh(αu)]/2$⁠, u is a scaled transmembrane voltage and v is a gating variable. The addition of weak diffusion in the gating variable (as well as replacing the Heaviside step function $Θ(u)$ with its smoothed version $Θα(u)$ and the additional term $Θα(v−1)(v−1)$⁠) serves to smooth out the unphysical discontinuities and singularities of the original model without noticeably changing its dynamics. This is done in order to facilitate computation of unstable solutions using the Newton-Krylov solver. It is well established, although perhaps not well known, that all relevant ions and small secondary messenger molecules pass through the gap junctions between cardiac cells.^33,34 Hence, inclusion of weak diffusion of the gating variable is justified physiologically.

Unless otherwise stated, the parameter values used were $u*=1.5415$⁠, M = 4, $ϵ=0.01, D11=1.1×10−3 cm2/ms,D22=D11/20$⁠, α  =  32, and $β=[1−exp (−R)]−1$⁠, where R = 1.273 is the restitution parameter. We also used no-flux boundary conditions $∇u·n̂=0$⁠, where $n̂$ 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 the pacing frequency is increased. Yet, with just two variables, as opposed to several tens of variables in the most complex models,³⁵ 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 L = 2.489 cm, which is close to the minimal size required for sustained SWC in the Karma model with the present choice of parameters. The PDE (1) was discretized in space using finite differences, with a five-point stencil for the Laplacian, on a 96 × 96 computational grid, which corresponds to the size of one computational cell $Δx=Δy=0.0262$ cm roughly equal to the physical length of one cardiac cell. The no-flux boundary conditions were approximated to second order in $Δx$ by the method of ghost points.³⁶ The values of u and v on this grid define the state of the system as a point in a 18432-dimensional state space. The discretized equations were solved using the 4th order Runge–Kutta method with time step $Δt=0.01$ ms chosen to be much smaller than the shortest characteristic time scale of the model, $τu=2.5$ ms.

The pacing was accomplished by applying a current to a 5 × 5 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 t_n the time at which the n-th pulse is applied, the pacing current can be written as

where $g(r)=1$ for paced grid points and zero otherwise.

We investigated stability of the excitation waves produced by pacing with gradually decreasing pacing interval $Tn=tn−tn−1, n=1,2,…$⁠, starting at $t=t0=0$⁠. 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).³⁷ However, for the Karma model considered here, neither of these conditions is necessary for wave breakup.

Let us define $T[i,j]$⁠, j > i, to be a sequence of pacing stimuli with a constant interval: $Ti=Ti+1=…=Tj$⁠. Consider the pacing protocol A: $T[1,10]=120$ ms, followed by $T[11,20]=110$ ms, $T[21,30]=100$ ms, $T[31,40]=90$ ms, and $T[41,80]=80$ ms (cf. Fig. 1) with the initial condition corresponding to the stable uniform equilibrium $u=[0,0]$⁠. The protocol B is identical to protocol A, except for $T[41,50]=85$ ms (dashed line in Fig. 1). For protocol B, each pacing stimulus generates a traveling wave with a 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 $2T=160$ 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, the 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 $t=t40$ 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 A. This is a signature of multistability often associated with subcritical instabilities; similar phenomenology is found during bypass transition to turbulence in fluid flows.^38,39

Let us start the analysis of the problem with a description of the various solutions to (1) for constant pacing, $T[1,∞]=T$⁠. In this case the dynamics can be reduced to a Poincaré map

that relates the states $un(r)=u(r,tn)$ at times t_n = nT (i.e., just before the nth pacing stimulus). In particular, time-periodic solutions $uk(r,t)$ of (1) with period kT, $k=1,2,…$⁠, correspond to k-cycles ${ukl(r)=uk(r,lT),l=1,…,k}$ of the Poincaré map, where $uk,l+1(r)=Φ[ukl(r)]$ for l < k and $uk1(r)=Φ[ukk(r)]$⁠. States $ukl(r)$ are fixed points of the k-th iterate $Φk$ of $Φ$⁠, and so they are solutions of the nonlinear equation

In particular, the 1:1 response $u1(r,t)$ corresponds to the fixed point of the Poincaré map $u11=Φ[u11]$⁠. At faster pacing rates, two different 2-cycles appear: ${u21,u22}$ corresponding to the alternans solution $u2(r,t)$ (2:2 response) and ${u′21,u′22}$⁠, corresponding to the conduction block solution $u′2(r,t)$ (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 $δu(r,t)=u(r,t)−uk(r,t)$ from the reference solution u_k. The dynamics of finite disturbances is described by

For infinitesimal disturbances, (7) reduces to a linear PDE

where $L=D∇2+J$ and $J(t)=∂uf|uk(t)$ 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 t_i to t_f > t_i in the tangent space,

can be found by integrating (8). In particular, the linear stability of $uk(r,t)$ is determined by the eigenvalues (Floquet multipliers) λ_ki of $U(kT,0)=∂Φk[u]/∂u$⁠,

where $eki(r)$ are the corresponding eigenfunctions (Floquet modes). In particular, if $|λki|<1$ for all i, any infinitesimal disturbance $δu(r,t)$ will eventually decay to zero, so that $u(r,t)→uk(r,t)$ for $t→∞$⁠. In the following discussion, we will index the eigenvalues in the order of decreasing absolute value, $|λk1|≥|λk2|≥|λk3|≥⋯$⁠.

Asymptotic stability, however, does not imply that initial disturbances will decay monotonically. Let us define the L₂ norm

where the scalar product is given by

and the integral is over the entire computational domain $0<x,y<L$⁠. The adjoint of the evolution operator with respect to the scalar product (12) defines the evolution of disturbances $δv$ in the adjoint of the tangent space,

backward in time (see, e.g., Ref. 40 for discussion and derivation of the adjoint equations). In practice, the action of $U†$ on a vector is most conveniently computed by solving the linear PDE

with “initial” condition defined at the final time t_f.

If the evolution operator U is nonnormal (i.e., $UU†≠U†U$⁠), some disturbances may grow transiently, before asymptotic decay takes over.²⁵ In the present case, the Jacobian J (and hence the differential operator $L$⁠) is non-self-adjoint due to coupling between the slow variable v and the fast variable u, which generates nonzero off-diagonal elements J₁₂ and J₂₁. As a result, the evolution operator U is indeed nonnormal, requiring the use of generalized linear stability theory^41,42 for describing transient response to perturbations associated with changes in the pacing rate.

Nonnormality and transient growth, however, are not unique to the model and geometry considered here. As we have shown previously using the Noble model⁴³ and the 4-variable Fenton model⁴⁴ in the context of one-dimensional paced cardiac tissue (Purkinje fibers), significant transient growth may arise, for example, due to the action of closed-loop feedback control aimed at stabilizing an otherwise unstable 1:1 solution. However, nonnormality (and hence transient dynamics) characterizes every known model of cardiac tissue, in any geometry, even without pacing or feedback. Some examples where non-self-adjointedness of the linearized evolution equations has been explicitly noted include the studies of spiral waves in the complex Ginzburg–Landau equation,⁴⁵ Fenton–Karma,⁴⁶ FitzHugh–Nagumo,⁴⁷ Barkley,⁴⁸ and Karma⁴⁹ models as well as scroll waves in generic excitable media^50,51 and the Beeler–Reuter–Pumir model of cardiac tissue.⁵² 

Given that time-periodic states $uk(r,t)$ could be stable or unstable, we computed them as solutions $uk1$ of (6) using a Newton–Krylov solver.^32,43 The spectrum of the evolution operator $U(kT,0)$ was computed using the Arnoldi method implemented by the Matlab function eigs⁵³ with accuracy yielding eight significant digits. In both instances, matrix-free evaluation of the matrix-vector product $U(t,0) δu$ was accomplished via simultaneous numerical time-integration of (1) and (8). The matrix-free calculation of the product $U(t,0)†δv$ involves backward integration and cannot be accomplished by integrating (1) and (14) simultaneously. Instead a precomputed and interpolated solution $uk(r,t)$ was used to compute the solution to (14) using the procedure described in Ref. 49. The time integrators were implemented as OpenCL kernels⁵⁴ 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.

The 1:1 solution $u1(r,t)$ exists for all T considered in this study, i.e., T > 50 ms. It is stable for $T>T2:2$ and becomes unstable for $T<T2:2$⁠, where $T2:2≈87.4$ ms for R = 1.273 (cf. Fig. 4). 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 $u2(r,t)$ with period 2 T is created. The bifurcation is supercritical (cf. Fig. 4), in agreement with the theoretical predictions,¹³ and the 2:2 solution exists for all $T<T2:2$ (at least down to T = 50 ms).

The 2:2 solution represents discordant alternans with four stationary nodal lines, as Fig. 5 illustrates. The alternans is created discordant, since the domain is sufficiently large, according to theoretical analysis.^13,14 In experiments^11,55 and simulations^29,37,56 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 the transition to fibrillation.

The stability balloon for the 2:2 solution in the (T, R) plane is shown in Fig. 6. This solution itself undergoes a subcritical Hopf bifurcation when the neutral stability curve $T=TH(R)$ (shown as a solid line) is crossed. The minimum of the neutral stability curve is at $Rc=1.279$ and T= 66.5 ms. The value R = 1.273 considered here is below R_c, so the 2:2 solution turns out to be linearly stable for all T. The corresponding Floquet spectra for different values of T are shown in Fig. 7. The change in the leading eigenvalue λ₂₁ associated with a decrease in T is shown in Fig. 7(a). Although this eigenvalue never leaves the unit circle, it approaches this circle for $T≈66.5$ 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. 6, as illustrated by Fig. 1 for example, 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 R = 1.273 and T = 80 ms. In addition, there is one more stable state $u′2(r,t)$ which describes a 2:1 response and exists for the same range of parameters. It is created in a saddle-node bifurcation as T decreases below a critical value $T2:1$ (⁠$T2:1=116.2$ ms for R = 1.273).

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 T = 80 ms that starts, respectively at t₄₀ or t₅₀). The perturbation $δuB0=uB(r,t50)−u21$ [cf. Figs. 8(c) and 8(d)] in protocol B eventually decays, in agreement with the prediction of linear stability analysis, while the perturbation $δuA0=uA(r,t40)−u21$ [cf. Figs. 8(a) and 8(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 T = 80 ms, the asymptotic regime is determined by the initial condition. In particular, there is a neighborhood $Ω2:2=Ω21∪Ω22$ of the 2-cycle ${u21,u22}$⁠, such that for all initial conditions $u0∈Ω2:2$⁠, the orbit $un$ approaches this 2-cycle as $n→∞$ (⁠$Ω2:2$ is known as the basin of attraction of the 2-cycle). Similarly, the 2-cycle ${u′21,u′22}$ has a basin of attraction $Ω2:1=Ω′21∪Ω′22$⁠. The basin $Ω2:2$ has a finite size in certain directions. In particular, the initial condition $uB0≡uB(r,t50)$ in protocol B lies inside $Ω2:2$⁠, so the asymptotic state corresponds to the 2:2 solution. On the other hand, the initial condition $uA0≡uA(r,t40)$ in protocol A lies outside both $Ω2:2$ and $Ω2:1$⁠, so the asymptotic state corresponds to SWC.

Since the evolution operator describing the linearized dynamics is nonnormal, the basin of attraction $Ω2:2$ becomes quite thin in some directions. Figure 9 shows a two-dimensional cartoon of this; the actual dimensionality of $Ω2:2$ 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 $||δu||$ from $u21$ to the boundary $∂Ω2:2$ 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 $∂Ω2:2$ that is the closest to the 2-cycle defines the critical optimal perturbation $δuoptcr$⁠. Although the shape of $Ω2:2$⁠, and hence the perturbation $δuoptcr$⁠, is defined by the nonlinear Poincaré map (5), as we illustrate below, both the direction of $δuoptcr$ and the strong transient growth associated with perturbations in this direction can be understood reasonably well within the linear approximation.

We will follow an analysis similar to that found in Refs. 25 and 57, but generalized for time-periodic solutions. Let us compute how strongly an initial disturbance $δu(r,0)$ is amplified as a result of transient growth. With the help of (9) and (11) we find

Since $U†(t,0)U(t,0)$ is self-adjoint, it has real eigenvalues $σi2(t)≥0$ and orthogonal eigenvectors $qi(r,t)$⁠,

with eigenvalues sorted from the largest to the smallest. Expanding the initial condition in the eigenvector basis yields

where $ηi(t)=〈qi(r,t),δu(r,0)〉$⁠, such that

The maximum of $||δu(r,t)||$ over all initial disturbances with a fixed norm $||δu(r,0)||$ is reached when all $ηi(t)=0$⁠, except for i = 1, i.e., $δuopt(r,0)∝q1(r,t)$⁠.

Let $||qi(r,t)||=1$ and define $pi(r,t)=U(t,0)qi(r,t)/σi(t)$⁠. The action of the evolution operator $U(t,0)$ can now be expressed in a simple manner. Multiplying both sides of (17) by $U(t,0)$ and recognizing that $δu(r,0)$ is arbitrary, we obtain

which is the singular value decomposition (SVD) of $U(t,0)$ with $σi(t), pi(r,t)$⁠, and $qi(r,t)$ being, respectively, the singular values, left singular vectors, and right singular vectors of $U(t,0)$⁠. The leading singular vectors and singular value have a useful interpretation: as we have already found, the right singular vector $q1(r,t)$ determines the shape of the optimal initial perturbation $δuopt(r,0)$ that is amplified the most [in the sense of the L₂ norm (15)] at time t. The left singular vector $p1(r,t)$ determines the shape this optimal perturbation acquires at time t. Finally, the singular value determines the magnitude of transient amplification $σi(t)=||δuopt(r,t)||/||δuopt(r,0)||$ which corresponds to this particular initial condition.

For small initial disturbances, the time at which the maximal transient amplification $σ1(t)$ 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 $maxtσ1(t)$ is numerically rather expensive, since SVD requires alternate direction integration of the linear PDEs (8) and (14) on the time interval $[0,t]$ and there is no simple relation that allows one to compute the singular values at time t + dt based on their values at time t. We therefore performed this calculation over many pacing periods on a coarse grid with $Δt=T/10$⁠.

The leading singular value $σ1(t)$ is shown in Fig. 10. 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, $σ1(t)$ appears to be a discontinuous function of time. In fact, this is not the case, as the time-resolved calculation over the interval $[26T,30T]$ illustrates (cf. Fig. 11). We find that $σ1(t)$ has a rather fine structure on a time scale smaller than one pacing period. The maxima of transient amplification have a pattern that is 2T-periodic; they are located at $t≈(2n+0.6)T$⁠, where n is an integer. The peaks' values, on the other hand, are modulated with a period $2π/arg(λ21)≈18.9T$ determined by the leading Floquet multiplier. The largest transient amplification (⁠$σ1=804$⁠) is achieved at $tmax=28.6T$ and the second highest value is achieved at $t=26.6T$⁠. Recall that for protocol A it takes a time interval of roughly 15 T until the dynamics transition from the neighborhood of the 2:2 state to SWC. These times are all of the same order of magnitude, suggesting that the memory effect extends over times scales of order $20T=1600$ ms, much longer than the naive estimate $ϵ−1τu=250$ ms would suggest.

The right singular vectors [cf. Figs. 12(a) and 12(b)] associated with the time instances when $σ1(t)$ is near the peak value $σ1(tmax)$ are indistinguishable, implying that their spatial structure, which determines the optimal initial perturbation $δuopt(r,0)$ is a robust feature of the dynamics. Indeed, the time-dependence of $σ1(t)$ for $t≳20T$ can be reproduced with extremely high accuracy by the norm of the disturbance which corresponds to the optimal one, e.g., $δuopt(r,0)=sq1(r,tmax)$⁠, where s is the amplitude of the initial disturbance. As Fig. 11 illustrates, the norm $||δuopt(r,t)||$ is indistinguishable from $σ1(t)$ once both have been normalized by their maximal values. Effectively, this means that the initial disturbance $sq1(r,tmax)$ is in fact optimal for all times $t≳20T$⁠. This optimal disturbance is localized around the pacing site (this is analogous to the localization of adjoint eigenfunctions of spiral wave solutions to the core region⁴⁹).

The corresponding left singular vector $p1(r,tmax)$ is shown in Figs. 12(c) and 12(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 L₂ norm (11), which involves integration over the entire spatial domain. Furthermore, the right singular vector is dominated by the v component, so that the dynamics are more sensitive to the perturbation of the v variable than the u variable. For the left singular vector the opposite is true, hence the perturbation affects the u variable to a much larger degree than the v 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.⁴⁹ 

As expected, generalized linear stability analysis correctly predicts the dynamics produced by optimal initial disturbances of sufficiently small amplitude s. Evolving the optimal initial perturbation with magnitude $s=10−6$ using the nonlinear Eq. (7), we find that $δuopt(r,t)$ grows transiently, achieving the predicted amplification $σ1(tmax)$ at time $tmax$ [cf. Fig. 13(a)]. For $t≳30T$ the asymptotic exponential decay $||δuopt(r,t)||∝eκt$ with the predicted rate $κ=ln|λ21|/(2T)$ sets in, where $|λ21|=0.8669$ is the leading Floquet multiplier. A stronger perturbation with $s=10−3$ 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.⁵⁸ In that case as well, transient growth of small localized perturbations was traced to the nonnormality of the evolution operator.⁵⁹ A 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.^43,44 In fact, increasing transient amplification was identified as the reason for breakdown of feedback control for long fibers controlled at the pacing site.

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 $u21$ that causes a transition to a different solution (⁠$u22, u′21, u′22$⁠, or SWC) is given by the distance from $u21$ to the nearest boundary of the corresponding basin of attraction (cf. Fig. 9). For linearly optimal initial perturbations, $δuopt(r,0)=sq1(r,tmax)$⁠, we find that the asymptotic state is $u21$ for $s≤1.227$⁠, as Fig. 13(b) shows. For $1.228≤s≤1.5949155$ the asymptotic state is $u22$⁠, which is also a 2:2 state, but with a different phase. For s = 1.5949156 the asymptotic state is again $u21$⁠, but for s = 1.5949157 the asymptotic state is SWC, which persists for at least 400 T. Increasing s further we find either persistent SWC or transient SWC that eventually gives way to one of the 2:1 states (⁠$u′21$ or $u′22$⁠).

The extreme sensitivity of the outcome to the choice of initial conditions is a characteristic feature of chaotic attractors that have fractal basin boundaries.⁶⁰ Since $u21$ and $u22$ are parts of the same 2-cycle, both $Ω21$ and $Ω22$ share a part of their boundary with the basin of attraction Ω_swc 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.⁵⁷ More recently, edge states have been studied extensively as gateways mediating subcritical transition to turbulence in shear fluid flows.^61,62

The 1:1 state $u11$ is an example of such an edge state. It has one unstable direction and its stable manifold forms the boundary between $Ω21$ and $Ω22$⁠. As Fig. 13(b) illustrates, the solutions corresponding to perturbations with s = 1.227 and s = 1.228 approach and closely follow the same periodic solution for $40T≲t≲160T$⁠, before eventually separating and approaching, respectively, $u21$ and $u22$⁠. This periodic solution is indeed $u11$⁠, as shown in Fig. 14. As expected, the two trajectories approach $u1(r,t)$ along the least stable direction and separate along the unstable direction, with the rates predicted by the linearization.

The basin boundary $∂Ωswc$ of the chaotic attractor lies further away from $u21$ than the basin boundary $∂Ω22$ of the 2:2 state $u22$⁠. An upper bound for the distance to $∂Ωswc$ is given by the norm $scr=||δuoptcr(r,0)||=1.5949157$ of the critical linearly optimal initial disturbance. While the generalized linear stability analysis may capture the qualitative features of the smallest initial disturbance required to trigger a transition from the 2:2 state to SWC, a fully nonlinear calculation is required to obtain a quantitatively accurate prediction for both its magnitude and shape.

Figure 15 shows the evolution of the initial condition $u(r,0)=u21+δuopt(r,0)$⁠, which corresponds to the linearly optimal initial disturbance with a magnitude s = 1.59492 that places it just outside of $Ω2:2$⁠. 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 $t=tr≈19.05T$⁠. Recall that, according to the 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 $U(tmax,0)$⁠: $δuopt(r,tmax)=σ1(tmax)sp1(r,tmax)$⁠. The location of conduction block is indeed found to be consistent with the shape of the left singular vector $p1(r,tmax)$⁠, which has maximal amplitude in the same corner of the domain, opposite the pacing site [cf. Figs. 12(c) and 12(d)]. However, the overall shape of the actual deviation $δuopt(r,tr)=u(r,t)−u2(r,t)$ (cf. Fig. 16) is noticeably different from $p1(r,tmax)$ since nonlinear terms become important for finite amplitude disturbances

While the first instance of conduction block generates a pair of spiral waves at $t≈19.35T$⁠, 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 $t≈20T$⁠. Another instance of conduction block occurs at $t≈22.85T$⁠, again generating two spiral waves at $t≈23.25T$⁠. 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 furthest 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 discussed previously, from $t=t40$ on, protocol A corresponds to constant pacing with T = 80 ms. The initial condition $uA0=u(r,t40)$ for the last stage of protocol A corresponds to a fairly significant perturbation away from the 2:2 state $u21$⁠. The norm $||δuA0||=||uA0−u21||=104.4$ of this perturbation is two orders of magnitude larger than $scr=||δuoptcr||$ (and is close to the norm of $u21$ itself). On the other hand, the spatial structure of this perturbation [cf. Figs. 8(a) and 8(b)] is quite different from that of the optimal perturbation [cf. Figs. 12(a) and 12(b)] and therefore it does not experience transient amplification, as Fig. 17(a) illustrates.

This however does not mean that transient dynamics do not play a role in the evolution of this perturbation; rather they manifest in a more subtle way. The component of $δuA0$ along $q1$ is 36 times smaller than the component along $p1$⁠, so transient growth of the perturbation along $q1$ is masked by its much larger orthogonal complement. As Fig. 17(b) shows, the magnitude of the projection $〈q1,δu〉$ stays at or below its initial value for around 30 T, after which it starts growing quickly, indicating the transition to SWC (recall that maximal transient amplification is predicted to occur at $28.6T$⁠). This result is consistent with the analysis we have performed for the critical optimal disturbance: while $|〈q1,δuA0〉|≈1.2$ is small relative to $|〈p1,δuA0〉|$⁠, it is comparable to the critical value s_cr. Notice that the evolution of $δuA0$ along $q1$ follows very closely that of the critical optimal disturbance; it takes nearly the same number of pacing intervals before the transition to SWC happens in both cases. This illustrates that the transition to SWC is controlled mainly by the projection of the initial disturbance on the optimal direction and not its L₂ norm.

Like $δuA0$⁠, the perturbation $δuB0$ that corresponds to protocol B has a spatial structure [cf. Figs. 8(c) and 8(d)] that is quite different from that of the optimal one and so it does not experience significant transient amplification either, as Fig. 17(a) illustrates. In this case, the component of $δuB0$ along $q1$ is 150 times smaller than the component along $p1$ or, in absolute terms, $|〈q1,δuB0〉|≈0.32$⁠, which is notably smaller than the critical value s_cr. This difference with protocol A turns out to be crucial: as Fig. 17(b) illustrates, the projection $〈q1,δu〉$ decays nearly exponentially (aside from some modulation) for protocol B, as the system returns to the discordant alternans state. Comparing the dynamics in the two protocols, we find that, on the qualitative level, it is the magnitude of the projection of the perturbation on the optimal direction that mainly determines whether that perturbation is inside or outside of the basin of attraction of the 2:2 state (cf. Fig. 9). Of course, a fully nonlinear analysis is required to make quantitatively accurate predictions.

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.^11,12,29 Theoretical analysis of alternans in one dimension predicts^13,18 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^12,63,64 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 furthest 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 the conduction block.

However, in a marked departure from the prevailing paradigm, the transition from discordant alternans to the 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 transient dynamics of disturbances away from the stable alternans state, which are associated with the nonnormality of the evolution operator. This nonnormality is a generic feature of all cardiac tissue models of the form (1), so our results should be quite representative.

Nonnormality requires the use of generalized linear stability analysis, which predicts, for the model studied here, that the optimal (most dangerous) disturbances correspond to perturbations of the gating variable near the pacing site. At the pacing interval of 80 ms, infinitesimal disturbances of the optimal shape were found to be transiently amplified by almost three orders of magnitude in terms of the L₂ norm. 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,⁶⁵ where transient amplification increases with the Reynolds number, which plays the role analogous to the pacing rate.

Considering the L₂ norm of a disturbance allowed us to detect and quantify transient dynamics and, in particular, identify the long temporal scales over which transients evolve. However, for typical scenarios involving various pacing protocols, it is not the L₂ norm of a disturbance that ultimately determines its fate; rather it is the projection of the disturbance on the optimal direction or, more specifically, the magnitude of the perturbation in the (slow!) gating variable near the pacing site. For a protocol with a decreasing pacing interval, the perturbation in the gating variable is controlled by the rate at which the pacing interval is decreased. For instance, in protocol B this perturbation is almost three times smaller than in protocol A, which leads to completely different outcomes. Generally, transition from alternans to SWC should happen sooner (at a longer pacing interval) if the pacing interval is decreased faster.

Although the transition from discordant alternans to SWC involves intrinsically nonlinear effects, such as the conduction block, we found that generalized linear stability analysis is capable of describing transient dynamics and some aspects of the transition semi-quantitatively. In particular, this analysis successfully predicts the spatial structure of the optimal disturbances, the magnitude of the critical disturbance required to provoke transition, as well as the time scales on which reentrant waves appear. Of course, it may be possible that transition to SWC could be triggered by nonlinearly optimal perturbations that are even weaker. 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.⁶⁶ 

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 of anisotropic tissue that are required to correctly describe intra- and extra-cellular potentials and in more complex geometries. 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) affects the dynamical description of transition from discordant alternans to fibrillation.

This material was 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.