ESAIM: M2AN 44 (2010) 1225–1238 DOI: 10.1051/m2an/2010028 ESAIM: Mathematical Modelling and Numerical Analysis www.esaim-m2an.org BIFURCATIONS IN A MODULATION EQUATION FOR ALTERNANS IN A CARDIAC FIBER Shu Dai 1, 2 and David G. Schaeffer 2 Abstract. While alternans in a single cardiac cell appears through a simple period-doubling bifurcation, in extended tissue the exact nature of the bifurcation is unclear. In particular, the phase of alternans can exhibit wave-like spatial dependence, either stationary or travelling, which is known as discordant alternans. We study these phenomena in simple cardiac models through a modulation equation proposed by Echebarria-Karma. As shown in our previous paper, the zero solution of their equation may lose stability, as the pacing rate is increased, through either a Hopf or steady-state bifurcation. Which bifurcation occurs first depends on parameters in the equation, and for one critical case both modes bifurcate together at a degenerate (codimension 2) bifurcation. For parameters close to the degenerate case, we investigate the competition between modes, both numerically and analytically. We find that at sufficiently rapid pacing (but assuming a 1:1 response is maintained), steady patterns always emerge as the only stable solution. However, in the parameter range where Hopf bifurcation occurs first, the evolution from periodic solution (just after the bifurcation) to the eventual standing wave solution occurs through an interesting series of secondary bifurcations. Mathematics Subject Classification. 35B32, 92C30. Received September 29, 2008. Revised July 26, 2009. Published online April 15, 2010. 1. Introduction An abnormal cardiac rhythm known as alternans, a bifurcation of the action potential duration (APD) of cardiac cells under rapid pacing stimuli, is believed to be one precursor of life-threatening ventricular fibrillation and sudden cardiac death [5,6,16]. The APD in a single paced cell has been modelled by the following restitution relationship [10,15]: (1.1) An+1 = f (DIn ), i.e. the APD only depends on the previous diastolic interval (DI), where An denotes the nth APD. In such a model alternans appears as a period-doubling bifurcation. In extended tissue composed of multiple cardiac cells, the action potentials generated by the stimuli will propagate in the tissue. Since the conduction velocity (CV) depends on the DI, the APD of the cell will also be a function of the cell’s position. In the case of one dimension, suppose we have a cardiac fiber of length ℓ, which is stimulated periodically at its x = 0 end, say with period B . Assume each stimulus successfully Keywords and phrases. Bifurcation, cardiac alternans, modulation equation. 1 2 Mathematical Biosciences Institute, The Ohio State University, Columbus, OH 43210, USA. sdai@mbi.osu.edu Department of Mathematics and Center for Nonlinear and Complex Systems, Duke University, Durham, NC 27708, USA. Article published by EDP Sciences c EDP Sciences, SMAI 2010  1226 S. DAI AND D.G. SCHAEFFER generates an action potential that propagates down the fiber. Let An (x) be the duration of the nth action potential at the position x along the fiber. If the pacing is rapid such that the basic cycle length B is below the critical value Bcrit , alternans are expected. To analyze position-dependent APD’s, in [4] Echebarria and Karma proposed the following as the basis for a multiscale expansion: An (x) = Acrit − δA + (−1)n an (x), (1.2) where Acrit is the APD when pacing with period B = Bcrit , δA is the average shortening of APD resulting from decreasing B below Bcrit and an (x) is the amplitude of alternans at the nth beat. It is assumed that an (x) varies slowly from beat to beat; thus one may regard it as the discrete values of a smooth function a(x, t) of smooth time t, i.e. an (x) = a(x, tn ) where tn = n · B for n = 0, 1, 2, . . ., the times when stimuli are applied. A weakly nonlinear modulation equation for a(x, t) was derived in [3] under approximation3, which after the nondimensionalization with respect to time, is given by  1 x 2 ∂t a = σa + ξ ∂xx a − w∂x a − a(x′ , t)dx′ − ga3 , (1.3) Λ 0 where σ is the bifurcation parameter, which is dimensionless and proportional to Bcrit − B; Λ, w, ξ are positive parameters, each having the units of length, that are derived from the equations of the cardiac model; and the nonlinear term −ga3 limits growth after the onset of linear instability. Neumann boundary conditions ∂x a(0, t) = 0, ∂x a(ℓ, t) = 0 (1.4) are imposed on (1.3). The trivial steady state solution a ≡ 0 of (1.3)–(1.4) loses its stability as σ increases. In the previous paper [2], we analyzed the eigenvalues of the linear operator that maps a function a(x) to  1 x 2 ′′ ′ a(x′ )dx′ , (1.5) ξ a − wa − Λ 0 subject to Neumann boundary conditions, and we concluded that the first bifurcation we observe as σ increases or Hopf bifurcation otherwise, where the critical value Λ−1 may be steady state bifurcation if Λ−1 < Λ−1 c c depends on ℓ, w and ξ. In this paper, extending our previous result, we simulate and analyze the bifurcation and hence competition of multiple modes appear. of (1.3) when Λ−1 is close to Λ−1 c The organization of this paper is the following: Section 2 reviews our previous result on the bifurcation from a simple eigenvalue. Section 3 shows the simulated results for the dynamics of the solution for Λ−1 around Λ−1 c . In Sections 4–5 we show that the dynamics of (1.3) may be derived from a reduced three dimensional flow. In Section 6 we perform a bifurcation analysis and relate it to the simulations. Incidentally, in a future paper we will show that for Λ−1 in a broader neighborhood of Λ−1 c , more complicated dynamics is obtained, possibly including chaos. 2. Linear stability: bifurcation from a simple eigenvalue We follow the nondimensionalization steps as in [2]: if we define dimensionless variables Λ̄ = Λ · w3 ξ −4 , x̄ = x · w ξ −2 , ℓ̄ = ℓ · w ξ −2 , then the operator in (1.5) can be written as w2 ξ −2 · La, where  x̄ d2 a da −1 La = − a(x̄′ )dx̄′ . − Λ̄ dx̄2 dx̄ 0 (2.1) (2.2) 3Equation (1.3) is rigorously valid only if the wave length of the solution a(x, t), which can be either a standing wave or travelling wave, is much larger than ξ. See [4] for more detail. 1227 ( real part of ) the eigenvalues BIFURCATIONS IN CARDIAC ALTERNANS 0 Re Ω 1 = Re Ω 2 Ω1 Ω0 Ω2 −2 Ω3 −4 Ω4 −6 −1 Λ 0 1 2 3 4 −1 Λc Figure 1. The evolution of the real parts of the first five eigenvalues Ω0 , Ω1 , . . . , Ω4 of the linear operator in (2.2) vs. Λ−1 , assuming ℓ = 6. All eigenvalues are real for Λ−1 sufficiently small, but, except for Ω0 , they become complex as Λ−1 increases. Λ−1 c ≈ 2.837 is the crossover the eigenvalue which has largest real part, Ωmax , is complex. point such that if Λ−1 > Λ−1 c And we rescale the time t and parameters σ and g t̄ = t · w2 ξ −2 , σ̄ = σ · w−2 ξ 2 , ḡ = g · w−2 ξ 2 . (2.3) In this notation, (1.3) may be rewritten as4 ∂t̄ a = σ̄a + La − ḡa3 . (2.4) For convenience below we will omit all the bars in (2.2) and (2.4). The boundary conditions are a′ (0) = 0, a′ (ℓ) = 0, (2.5) where ℓ is the dimensionless length of the cardiac fiber. Obviously a(x, t) ≡ 0 is a trivial solution to (1.3). When σ is small (including all negative values), the zero solution is linearly stable. As σ increases to beyond some threshold, the zero solution loses its stability. To investigate the bifurcation, we consider the eigenvalues of the linear operator (2.2). Suppose Ω0 , Ω1 , Ω2 , . . . are those eigenvalues and let Ωmax be the one having the largest real part. If σ + Re Ωmax < 0, the zero solution is linearly stable; otherwise it is unstable. This bifurcation at σ = −Re Ωmax can be either steady state or Hopf, depending on whether Ωmax is real or complex. The dimensionless linear operator (2.2), together with the boundary conditions (2.5) has two parameters ℓ and Λ. In our previous paper [2], we found that for ℓ sufficiently large, as Λ−1 varies, the eigenvalues of (2.2) behave as illustrated in Figure 1. In the figure ℓ = 6 and only the real parts of all eigenvalues are graphed (Remark: In [2] we computed the eigenvalues asymptotically for large ℓ). By the critical value of Λ−1 , say Λ−1 c , as denoted in the graph, we mean the crossover point where the eigenvalue with largest real part Ωmax changes from real to complex. We observe that all eigenvalues lie in the left half of the complex plane, i.e. they all have negative real parts. When Λ−1 < Λ−1 c , Ωmax = Ω0 . Since Ω0 is real, as σ increases from 0, we will first encounter a steady state bifurcation at σ = |Re Ωmax | = −Ω0 , and as we compute below beyond the bifurcation point the solution to equation (1.3) has a stable pattern of standing waves. On the other hand if Λ−1 > Λ−1 c , we have Ωmax = Ω1,2 , which is a complex pair and as σ increases from 0 we will encounter a Hopf bifurcation 4We could remove the factor ḡ by rescaling on the amplitude a(x, t), but little would be gained by this. In this paper, ḡ is fixed at 40 if not specifically mentioned. 1228 S. DAI AND D.G. SCHAEFFER 3.8 standing wave 3.6 (3) c 3.4 (2) σ d 3.2 b (1) 3 a 2.8 zero steady state (0) 2.6 2.6 2.7 2.8 Λ−1 c 2.9 3 3.1 3.2 3.3 3.4 −1 Λ Figure 2. Regions in parameter space where the solution of (1.3) exhibits different longterm behavior, for Λ−1 near Λ−1 when ℓ = 6. There are four regions with the following c behaviors: (0) trivial zero steady state solution, (1) pure periodic solution, (2) mixed-mode periodic solution and (3) standing wave solution. We also mark four points (a–d) in the graph, which correspond to the simulations in Figure 3 (a–d). The vertical arrows will be useful for the discussion below. at σ = |Re Ωmax | = −Re Ω1 , and at least for a small range of σ beyond the bifurcation point the solution has a pattern of travelling wave5. A calculation in [2] shows that for sufficiently large ℓ, to the leading order, the dimensionless critical value equals √ 71 + 17 17 −1 + O(ℓ−2 ). (2.6) Λc = 64 In [2] we discussed two specific models (see Tab. 5.1 in [2]), which show Λ−1 can be either larger or smaller −1 ≈ 48.4, while for a cable of the length in this than Λ−1 c . For the two-current model [13], the dimensionless Λ example, the dimensionless Λ−1 ≈ 2.67. As expected a travelling pattern occurs. For the Noble model [14], c ≈ 2.27. As expected a pattern of standing wave emerges. Values of Λ−1 between these Λ−1 ≈ 0.23 while Λ−1 c , may be anticipated from other ionic models. extremes, including values close to Λ−1 c In the critical case, i.e. when Λ−1 = Λ−1 , we have Ω = Re Ω , both modes become unstable simultaneously. 0 1,2 c In the next section we investigate numerically the behavior of the system near such a degenerate point. 3. Competition between modes near the critical point In the following sections we fix the dimensionless length of the cardiac fiber to be ℓ = 6 to illustrate the dynamics. We find by computation that for ℓ = 6, the critical point Λ−1 c ≈ 2.837. Figure 2 illustrates simulation results for the modulation equation (1.3) when Λ−1 is near the critical value Λ−1 c . Figure 3 shows four typical simulated solutions whose parameters correspond to points a-d in Figure 2, which belong to regions (0–3) respectively. If (Λ−1 , σ) lies in the quasitriangular region (3) in Figure 2, as in Case d, the stable solution of (1.3) is a steady pattern. How this steady pattern appears as σ is increased (quasi-statically) depends on the value of Λ−1 (assumed fixed). For Λ−1 < Λ−1 c , as indicated by the vertical arrow on the left in Figure 2, this steady solution appears through a single, simple bifurcation. However 5The instability can either be convective or absolute. In the limit ℓ → ∞, the spectrum {Ω } for n ≥ 1 is nearly continuous. n We have convective instability if the group velocity is nonzero, and absolute instability otherwise. See [4] for a more complete discussion. 1229 BIFURCATIONS IN CARDIAC ALTERNANS 0.2 0.1 0.1 a(x2 ,t) a(x,t) (a) Λ−1 = 3.3, σ = 2.9 0.2 0 −0.1 −0.2 0 −0.1 0 1 2 3 −1 (b) Λ 4 5 6 −0.2 −0.2 −0.1 0.2 0.2 0.1 0.1 0 0 −0.1 −0.1 −0.2 0 1 2 3 −1 (c) Λ 4 5 6 −0.2 −0.2 −0.1 0.2 0.1 0.1 0 0 −0.1 −0.1 0 1 2 3 −1 (d) Λ 4 5 6 −0.2 −0.2 −0.1 0.2 0.1 0.1 0 0 −0.1 −0.1 0 1 0.2 0 0.1 0.2 0 0.1 0.2 0 0.1 0.2 = 3.3, σ = 3.6 0.2 −0.2 0.1 = 3.3, σ = 3.45 0.2 −0.2 0 a(x1 ,t) = 3.3, σ = 3.2 2 3 4 5 6 −0.2 −0.2 −0.1 position on the fiber Figure 3. Four typical simulations of the amplitude equation (2.4), with g = 40 fixed and other parameters corresponding to the points a−d in Figure 2, and the initial condition is given by a(x, 0) ≡ 0.1. The left column shows the graph of a(x, t) for t = 36, 37 and 38 (solid, dashed and dash-dot). For case (a), the zero steady state, and case (d), the nonzero standing wave, the difference among the three curves is negligible. The right column is the phase plane portrait, after the transient, obtained by plotting a(x1 , t) versus a(x2 , t), where x1 = 4.5 and x2 = 5.1 are two fixed positions on the fiber. if Λ−1 > Λ−1 c , as indicated by the arrow on the right in Figure 2, the evolution to a steady pattern proceeds through several intermediate states. The first bifurcation is to a pure periodic solution (Case b), followed by a secondary bifurcation to a mixed periodic/stationary solution (Case c), and finally bifurcation to the steady solution. Let us discuss further the two nonsteady cases (b and c). In both cases the solution has a travelling pattern, as shown in the left column of Figure 3. Graphs in the right column in the figure clarify the difference between Cases b and c. In this column fixed points x1 = 4.5 and x2 = 5.1 are chosen and a(x1 , t) is plotted versus a(x2 , t). In Case b, the average of a(x, t) over time is nearly zero, while in Case c this average differs significantly from zero. To better understand the above behavior, consider the solution a(x, t) at a fixed point x1 = 4.5 as σ increases and Λ−1 = 3.3, i.e., following the arrow on the right in Figure 2. Figure 4 shows how the quantity max a(x1 , t) − min a(x1 , t) t t (3.1) 1230 S. DAI AND D.G. SCHAEFFER 0.3 0.2 0.1 (0) (1) 0 −0.1 (3) (2) Max − Min Average 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 σ Figure 4. The simulation result for the oscillation amplitude, maxt a(x1 , t) − mint a(x1 , t) (dashed curve), and the average of a(x1 , t) (solid curve) in one period for various values of σ when Λ−1 = 3.3 and x1 = 4.5. The four labels refer to regions in Figure 2: (0) for σ < 3.14, zero steady state solution, both amplitude and average are zero, (1) for 3.14 < σ < 3.32, pure periodic solution, the amplitude is nonzero but the average is zero (to leading order), (2) for 3.32 < σ < 3.49, mixed mode solution, both amplitude and average are nonzero, (3) for σ > 3.49, standing wave, with no oscillation and nonzero average. (dashed line) and the average  1 T a(x1 , t) dt, (3.2) T 0 where T is the period (solid line), vary with σ for 3 < σ < 3.7. In region (1), which we call the pure-periodic region, the average of the solution is close to zero. By contrast, in region (2), although the solution continues to oscillate, its average is nonzero. Finally, in region (3) the solution is steady (and equal to its average). In the following sections we provide supporting theory for the above observation. In particular we shall see that the transitions between region (1) and (2) and between (2) and (3) represent secondary bifurcations that appear from unfolding the degenerate bifurcation for Λ−1 = Λ−1 c . 4. Reduction to a three dimensional system We consider the case when the parameter Λ−1 exactly equals to the critical value Λ−1 c , so we have Ω0 = Re Ω1,2 > Re Ω3,4 > . . . , which are all negative. For appropriate σ, we apply center-manifold theory [1,9] to show that the dynamics of the solution to (1.3) is determined by the first three modes corresponding to Ω0,1,2 . First we rewrite the modulation equation (1.3) in a more convenient way. Define σ̃ = σ − σc , the increment of σ away from its bifurcation value σc = −Ω0 . In dimensionless variables, (1.3) can be rewritten as ∂t a = σ̃a + L a − a3 , where L a = σc a + a′′ − a′ − Λ−1 c  (4.1) x a(x′ )dx′ 0 and we have assumed the positive nonlinear coefficient g = 1 by a scaling a → g −1/2 · a. (4.2) BIFURCATIONS IN CARDIAC ALTERNANS 1231 Table 1. Values of the coefficients hji0 i1 i2 in system (4.14). (i0 , i1 , i2 ) 3, 0, 0 2, 1, 0 2, 0, 1 1, 2, 0 1, 1, 1 1, 0, 2 0, 3, 0 0, 2, 1 0, 1, 2 0, 0, 3 h0i0 i1 i2 0.168 −0.070 −0.346 0.081 0.010 0.582 0.004 −0.053 −0.021 −0.084 h1i0 i1 i2 0.008 0.132 −0.152 −0.104 0.021 0.160 0.107 −0.228 0.262 −0.040 h2i0 i1 i2 0.055 −0.024 −0.188 0.012 0.017 0.146 −0.015 0.023 0.012 0.096  n = Ωn + σc = Ωn − Ω0 for n = 0, 1, 2, . . . Thus the first three eigenvalues are all The eigenvalues of L are Ω on the imaginary axis:  0 = 0, Ω  1,2 = ±i ω, Ω (4.3) where ω > 0 is real. All other eigenvalues are in the left half plane. Let φ0 (x) be the eigenfunction with  1,2 , i.e.  0 and φ1 (x) ± i φ2 (x) be the eigenfunctions with eigenvalues Ω eigenvalue Ω L φ0 = 0 and L (φ1 ± i φ2 ) = ±i ω(φ1 ± i φ2 ), (4.4) where φ0 , φ1 and φ2 are all real. Let E c = span{φ0 , φ1 , φ2 }, i.e. E c is the central subspace. We study the dynamics of a(x, t) through its projection onto the central subspace. We regard the solution to (4.1) as a flow in the function space of L2 on the interval (0, ℓ), i.e. square integrable functions. By the central manifold theorem [1,9], there exists a central manifold M which is invariant under the flow and tangent to the central space E c , and there is also a neighborhood of zero in L2 , denoted by U , with the following property: if a(·, t) is a solution of (4.1) such that for all t > 0, a(·, t) ∈ U , then the distance from a(·, t) to M converges to zero exponentially. In other words, to understand the long-time dynamics of (4.1) near 0, it is sufficient to examine the flow on M. The flow on M is a three dimensional ODE, which can be formulated by using the coordinates of E c . We first introduce the adjoint operator of L , which is defined as below: L ∗ a = σc a + a′′ + a′ − Λ−1 c  ℓ a(x′ )dx′ , (4.5) x with boundary conditions a′ (0) + a(0) = a′ (ℓ) + a(ℓ) = 0. The adjoint operator L of L ∗ , i.e., ∗ (4.6) has the same eigenvalues as L . Let ψ0 , ψ1 ± i ψ2 be the first three eigenfunctions L ∗ ψ0 = 0 and L ∗ (ψ1 ± i ψ2 ) = ∓i ω · (ψ1 ± i ψ2 ), (4.7) where ψ0 , ψ1 and ψ2 are all real.We may impose the following conditions of orthogonality: ψi , φj = δij , for i, j = 0, 1, 2, (4.8) where the inner product , is taken in the L2 -sense. Figure 5 shows a possible choice of φi (x) and ψi (x) for i = 0, 1, 2 by numerics. To obtain the eigenfunctions φi ’s, we discretized L to second order with space step dx = 0.003, where the derivatives are approximated by finite differences and the integral term is approximated by trapezoidal rule. Then we use LAPACK to find the eigenvalues of the resulting matrix, the first three of which are just 0, ±i ω as in (4.4). We solve the linear system L φ0 = 0 to find the eigenfunction φ0 . We also solve the complex system L (φ1 + i φ2 ) = i ω(φ1 + i φ2 ) by separating into its real and imaginary parts. The eigenfunctions φi ’s in Figure 5 are obtained after L2 -normalization: ℓ 2 ℓ 2 2 0 φ0 dx = 1 and 0 (φ1 + φ2 ) dx = 1. The eigenfunctions ψ0,1,2 of the adjoint operator are obtained in a similar way. However, they are not normalized in L2 , but by requiring orthogonality (4.8). 1232 S. DAI AND D.G. SCHAEFFER 1 0.5 φ1 15 φ2 5 0 ψ2 0 −5 −0.5 −1 ψ0 10 φ0 0 1 2 3 4 position on the fiber 5 ψ1 −10 6 −15 0 1 2 3 4 position on the fiber A 5 6 B Figure 5. Illustration of the eigenfunctions φ0,1,2 in (4.4) and ψ0,1,2 in (4.7), assuming the length of the fiber ℓ = 6. They satisfy the conditions of orthogonality (4.8). We now parameterize M by E c . Let H = {ψ0 , ψ1 , ψ2 }⊥ ; then H is a complement of E c in L2 , i.e. L2 ∼ = E c ⊕H. 2  Let π be the projection onto E c with kernel H. Since M is tangent to E c , for all u = ui φi in E c ∩ U , there is i=0 a unique J(u) ∈ M such that πJ(u) = u. Moreover, the difference R(u) = J(u) − u belongs to H and satisfies R(u) = O(|u|2 ). (4.9) Suppose a(x, t) is a solution of (4.1) such that for all time a(·, t) ∈ M ∩ U . Let u(t) = πa(·, t). Of course a(·, t) = J(u(t)). To derive an ODE for u(t), we calculate u̇ = π ȧ = σ̃u + πL J(u) − π[J(u)3 ] =σ u + πL u + πL R(u) − π(u3 + O(|u|4 )). (4.10) Clearly πL u = L u for L u ∈ E c since u ∈ E c . Since R(u) ∈ H, for j = 0, 1, 2, we have ψj , L R(u) = L ∗ ψj , R(u) = 0, i.e. L R(u) ∈ H and hence πL R(u) = 0. Thus (4.10) can be rewritten as u̇ = σ̃u + L u − π(u3 ) + O(|u|4 ). (4.11) Let us regard E c , a subspace of L2 , as a three-dimensional vector space with coordinates (u0 , u1 , u2 ) defined by u= 2  ui (t)φi (x). i=0 To expand (4.11) in coordinates, we take the inner product of this equation with each ψj for j = 0, 1, 2. Since u3 − π(u3 ) ∈ H = {ψ0 , ψ1 , ψ2 }⊥ , we conclude that ψj , π(u3 ) = ψj , u3 for j = 0, 1, 2. By (4.4) and (4.8), we obtain the ODE system for each coordinate ui (t): ⎧  2 3  ⎪ ⎪ ⎪ u̇ u φ = σ̃u − ψ , + O(|u|4 ), i i 0 0 0 ⎪ ⎪ i=0 ⎪ ⎪  ⎨ 2 3  u̇1 = σ̃u1 − ωu2 − ψ1 , + O(|u|4 ), ui φi ⎪ i=0 ⎪ ⎪  ⎪ 2 3 ⎪  ⎪ ⎪ u̇2 = σ̃u2 + ωu1 − ψ2 , ui φi + O(|u|4 ). ⎩ i=0 (4.12) BIFURCATIONS IN CARDIAC ALTERNANS We expand  2  i=0 ui φi 3 1233 in each equation in (4.12), and the coefficient for each term ui00 ui11 ui22 , where i0 + i1 + i2 = 3, in the j-th equation is given by hji0 i1 i2 = ψj , φi00 φi12 φi03 =  ℓ ψj φi00 φi12 φi03 dx. (4.13) 0 ℓ For instance, h12,1,0 = ψ1 , φ20 φ1 = 0 ψ1 φ20 φ1 dx. Choosing the functions φi ’s and ψj ’s as in Figure 5, we computed all the coefficients numerically, and these are in Table 1. Therefore the reduced system (4.12) can be rewritten as the following:  ⎧ u̇0 = σ̃u0 − h0i0 i1 i2 ui00 ui11 ui22 + h.o.t., ⎪ ⎪ ⎨  u̇1 = σ̃u1 − ωu2 − h1i0 i1 i2 ui00 ui11 ui22 + h.o.t., ⎪ ⎪  ⎩ u̇2 = σ̃u2 + ωu1 − h2i0 i1 i2 ui00 ui11 ui22 + h.o.t., (4.14) where the summations are over nonnegative integers i0 , i1 , i2 such that i0 + i1 + i2 = 3 and h.o.t. means higher order terms. We rewrite (4.14) in the compact form u̇ = σ̃u + ω · T u + H(u) + h.o.t., where ⎛ 0 Tu = ⎝ 0 0 and H(u) is the vector-valued homogeneous cubic (4.15) ⎞ ⎛ ⎞ 0 0 u0 0 −1 ⎠ · ⎝ u1 ⎠ (4.16) 1 0 u2 polynomials of u0 , u1 , u2 in (4.14), including the minus sign. 5. Derivation of the normal form 5.1. Elimination of the nonresonant terms To investigate the dynamics of the reduced system (4.14), following [9] we perform a polynomial transformation of coordinates to obtain its normal form. Let H3 be the space of homogeneous polynomials of degree 3. We − → can regard H(u) in equation (4.15) as an element in the space H3 ⊕ H3 ⊕ H3 = H3 , a basis of which is given by ⎧⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎫ 1 0 0 ⎨ ⎬  3 2  u1 , u1 u2 , u21 u0 , u1 u22 , u1 u2 u0 , u1 u20 , u32 , u22 u0 , u2 u20 , u30 ⊗ ⎝ 0 ⎠ , ⎝ 1 ⎠ , ⎝ 0 ⎠ . (5.1) ⎩ ⎭ 0 0 1 Consider a transformation of the form (u0 , u1 , u2 ) = (v0 , v1 , v2 ) + P (v0 , v1 , v2 ), (5.2) − → where P ∈ H3 is a vector-valued homogeneous cubic polynomial. We substitute (5.2) into (4.15) to find v̇ = σ̃v + ω · T v + H(v) + adT (P )(v) + h.o.t. . (5.3) − → − → Here the adjoint operator adT (·) : H3 → H3 in above equation is defined by − → ad T (P )(v) = T P (v) − (DP ) · T v, ∀P ∈ H3 , (5.4) 1234 S. DAI AND D.G. SCHAEFFER Table 2. Values of the coefficients in the computed normal form (5.9). a1 a2 b1 b2 c1 c2 0.083 −0.096 −0.130 −0.422 0.003 −0.090 where DP = (∂j Pi ) is the 3 × 3 matrix. We shall write H(v) = H1 (v) + H2 (v), where H1 ∈ Ker(ad T ) and H2 ∈ Range(ad T ). Then by an appropriate choice of P such that ad T (P ) = −H2 , the cubic terms H(v) + ad T (P )(v) will reduce to H1 (v), i.e. the projection of the H(v) onto the kernel Ker(ad T ). To carry out this reduction, we let ad T (·) defined in (5.4) act on each vector of the basis (5.1) and expand the result in the same basis to find ad T (·) has the following matrix ⎛ S10 ad T = ⎝ O O O S10 I10 ⎞ O −I10 ⎠, S10 (5.5) where I10 is 10 dimensional identity matrix, O is the 10 × 10 zero matrix and ⎛ S10 ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ =⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 3 0 0 0 0 0 0 0 0 −1 0 0 2 0 0 0 0 0 0 0 0 0 0 0 −2 0 0 0 0 −1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 ⎞ 0 0 0 0 0 0 0 0 ⎟ ⎟ 0 0 0 0 ⎟ ⎟ −3 0 0 0 ⎟ ⎟ 0 −2 0 0 ⎟ ⎟. 0 0 −1 0 ⎟ ⎟ 0 0 0 0 ⎟ ⎟ 0 0 0 0 ⎟ ⎟ 0 0 0 0 ⎠ 0 0 0 0 (5.6) Using Maple to investigate matrix (5.5), we find that ad T is diagonalizable (in the complex sense) and its null − → space is a six dimensional subspace of H3 spanned by the following eigenvectors: ⎞ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ 3 ⎞ ⎛ r2 v0 0 0 0 v0 0 ⎝ 0 ⎠ , ⎝ 0 ⎠ , ⎝ r2 v1 ⎠ , ⎝ v1 v02 ⎠ , ⎝ −r2 v2 ⎠ , ⎝ −v2 v02 ⎠, v1 v02 v2 v02 r2 v1 r2 v2 0 0 ⎛ (5.7)  where for convenience we defined r = r(v) = v12 + v22 . On the other hand, the remaining twenty-four eigenvectors of ad T corresponding to nonzero eigenvalues span the range of ad T and can therefore be transformed away by an appropriate choice of P in (5.2). Thus the system (4.12) can be transformed to the following simplified form: ⎧ ⎨ v˙0 = σ̃v0 + a1 r2 v0 + a2 v03 + h.o.t., v˙1 = σ̃v1 − ωv2 + b1 r2 v1 + b2 v1 v02 − c1 r2 v2 − c2 v2 v02 + h.o.t., (5.8) ⎩ v˙2 = σ̃v2 + ωv1 + b1 r2 v2 + b2 v2 v02 + c1 r2 v1 + c2 v1 v02 + h.o.t. Using Maple to perform the computations, we find that the coefficients in Table 1 lead to a reduced system (5.8) with the coefficients a1,2 , b1,2 and c1,2 given in Table 2. 1235 BIFURCATIONS IN CARDIAC ALTERNANS 5.2. Scaling of the remaining terms We introduce polar coordinates such that v1 = r cos θ and v2 = r sin θ and for convenience we define z = v0 . Then the system (5.8) can be rewritten as the following ⎧ ⎨ ż = z(σ̃ + a1 r2 + a2 z 2 ) + O(|r, z|4 ), ṙ = r(σ̃ + b1 r2 + b2 z 2 ) + O(|r, z|4 ), ⎩ θ̇ = ω + O(|r, z|2 ). (5.9) We study the reduced bifurcation problem in the variables r and z. By scaling these variables z̃ = z ·  √ −a2 , r̃ = r · −b1 , (5.10) we may reduce this subsystem of (5.9) to the normal form  d dt z̃ = z̃(σ̃ − mr̃2 − z̃ 2 ) + O(|r̃, z̃|4 ), d dt r̃ = r̃(σ̃ − nz̃ 2 − r̃2 ) + O(|r̃, z̃|4 ), where m= (5.11) b2 a1 ≈ −0.64, and n = ≈ 4.40. b1 a2 (5.12) Since to lowest order θ̇ is a constant, an equilibrium of (5.11) with r = 0 corresponds to a periodic solution of (5.9). 6. Analysis of the bifurcation In previous sections for the case Λ−1 = Λ−1 exactly, we derived equation (5.9) for the restriction of (4.1) c to the center manifold. Let us show that equation (5.9) undergoes both a Hopf bifurcation and steady state bifurcation as σ̃ passes through zero. Motivated by [8,12], we expect that interactions between the two modes can lead to secondary bifurcations. Although later we shall specialize to the values of m and n in (5.12), any values satisfying n > 1, mn < 1 yield qualitatively similar behavior (see Chap. 10, [7]). Neglecting the h.o.t. in (5.11) and dropping the tildes on r̃, z̃, we obtain the system  ż ṙ  = f (z, r, σ̃) =  z(σ̃ − mr2 − z 2 ) r(σ̃ − nz 2 − r2 )  . (6.1) Note that this system has Z2 ⊕ Z2 symmetry under the transformations z → −z and r → −r. This symmetry helps in the enumeration of equilibria f (r, z, σ̃) = 0 of (6.1). Specifically, we have the four possible cases: • the trivial zero solution (0, 0); √ • a pure z-mode (± √σ̃, 0), corresponding to a standing wave solution; • a pure r-mode (0, σ̃, corresponding to a periodic solution of (5.9); • mixed mode which satisfies (6.2) σ̃ = mr2 + z 2 = nz 2 + r2 . The pure z-mode and the pure r-mode appear only while σ̃ ≥ 0, i.e., both modes bifurcate supercritically. Regarding possible mixed modes, given m and n as in (5.12) there is no value of σ̃ for which (6.2), a pair of linear equation in r2 and z 2 , has real nonzero solutions. Checking the eigenvalues of df (z, r, σ̃), we find when σ̃ < 0, the trivial zero solution is stable; and when σ̃ > 0, only the pure z-mode is stable. We summarize this information in the schematic bifurcation diagram Figure 6, where a solid line means stable and dashed line means unstable. 1236 S. DAI AND D.G. SCHAEFFER r − mode (periodic) z − mode (steady state) 0 zero mode Figure 6. Schematic bifurcation diagram to the system (5.11), where the horizontal direction denotes the value of σ̃ . The z-mode represents a symmetric pair of steady solutions of (5.9), and the r-mode represents a periodic solution of (5.9). A solid line means the mode is stable, and a dashed line means it is unstable. More generally, we assume Λ−1 in (1.3) is close to, but not exactly equal to, Λ−1 c . As illustrated in Figure 1, changing Λ−1 shifts both Ω0 and Re Ω1,2 . To model these shifts we consider an unfolding [7,11] of (6.1) with two auxiliary parameters µ1 and µ2 :     z(σ̃ − µ1 − mr2 − z 2 ) ż . (6.3) = r(σ̃ − µ2 − nz 2 − r2 ) ṙ Letting σ̂ = σ̃ −µ2 and µ = µ1 −µ2 , we may eliminate an inessential parameter and rewrite (6.3) in the following form:     ż (σ̂ − µ)z − mr2 z − z 3 . (6.4) = F (z, r, σ̂, µ) = σ̂r − nrz 2 − r3 ṙ As above, we can enumerate equilibria of (6.4) – solutions of F (z, r, σ̂, µ) = 0 – in four cases: • the trivial zero solution (0, 0); √ • for σ̂ > µ, a pure z-mode (± √σ̂ − µ, 0), corresponding to a standing wave solution; • for σ̂ > 0, a pure r-mode (0, σ̂, corresponding to a periodic solution; • a mixed mode which satisfies σ̂ = mr2 + z 2 + µ = nz 2 + r2 . The stability of the bifurcating solutions and the existence of (real) mixed-mode solutions depend on the sign of µ. Case I: µ < 0. In this case there are no mixed-mode solutions for any value of σ̂. Thus the equilibria of (6.4) trace out a bifurcation diagram as sketched in Figure 7A. Exchange of stability suggests that the trivial solution for σ̂ < µ and the z-mode for σ̂ > µ are stable and other solutions are unstable. This is easily confirmed from the Jacobian of (6.4). Case II: µ > 0. In this case, the mixed-mode equations may be solved to obtain r2 = (n − 1)σ̂ − nµ , mn − 1 In the range δ1 < σ̂ < δ2 , where z2 = (m − 1)σ̂ + µ · mn − 1 (6.5) nµ µ and δ2 = , (6.6) 1−m n−1 both right-hand-sides of (6.5) are positive. Thus (real) mixed-mode solutions exist for this range of σ̂. Since z → 0 as σ̂ → δ1+ and r → 0 as σ̂ → δ2− , the mixed-mode solution branch connects the pure r- and z-modes, as sketched in Figure 7B. Exchange of stability suggests the stability assignments of the figure, and these are readily verified. Speaking teleologically, we may say that the secondary bifurcations of the mixed-mode solution are needed for the primary branches to have the same stability as in Figure 6 for large σ̂. These bifurcation diagrams agree with the simulations. For fixed Λ−1 < Λ−1 c , when we increase σ from 0, the simulations in Figure 2 show a steady state bifurcation from the trivial solution in region (0) to a standing wave solution in region (3), which is the behavior in the bifurcation diagram Figure 7A. For fixed Λ−1 > Λ−1 c , when we increase σ from 0, Figure 2 shows we first encounter a bifurcation from the trivial solution in region (0) δ1 = 1237 BIFURCATIONS IN CARDIAC ALTERNANS z − mode (steady state) µ< 0 1.5 mixed mode 1.5 r − mode (periodic) µ >0 1 1 0.5 0.5 0 −2 z − mode (steady state) r − mode (periodic) 0 −1.5 µ −1 −0.5 −1 0 0 0.5 1 1.5 −1 −0.5 −1 00 δ1 0.5 µ δ 21.5 1 2 −1 B: Λ−1 > Λ c A: Λ < Λ c Figure 7. Schematic bifurcation diagrams for two cases Λ−1 < Λ−1 and Λ−1 > Λ−1 c c , where the horizontal direction denotes the value of σ̂. The solid line means the mode is stable and the dashed line means it is unstable. to a pure r-mode (periodic solution with zero average, to leading order) in region (1), and then it bifurcates to a mixed mode (periodic solution with nonzero average) in region (2), finally it bifurcates to the pure z-mode (standing wave solution) in region (3), all of which can be observed from Figure 7B as we increases σ̂. The above data also provide a check on the accuracy of truncating all h.o.t. We can obtain δ1 and δ2 from the simulation in Figure 4, where Λ−1 = 3.3. δ1 represents the increase in σ̂, or equivalently of σ, from the first bifurcation (emergence of the pure r-mode) to the second bifurcation (the appearance of the mixed mode), and similarly δ2 represents the increase in σ̂ or σ from the first bifurcation to the third. We find δ1 sim ≈ 3.32 − 3.14 = 0.18, δ2 sim ≈ 3.49 − 3.14 = 0.35. (6.7) On the other hand, since µ = µ1 − µ2 , where to the leading order µ1 and µ2 are shifts of the threshold values of σ̃ for the emergence of pure z-mode and r-mode respectively, we obtain µ = Re Ω1 − Ω0 . By computation of the eigenvalues, as shown in Figure 1, when Λ−1 = 3.3, we have µ ≈ 0.28. Thus the formula in (6.6) gives the theoretical result: nµ µ ≈ 0.17, δ2 theo = ≈ 0.36 (6.8) δ1 theo = 1−m n−1 for the values of m and n in (5.12), which matches the simulated result (6.7) well. 7. Conclusion We have studied the bifurcations of the modulation equation (1.3) for APD alternans propagating on a cardiac fiber. We observed the solutions undergo either a steady state or a Hopf bifurcation; which occurs first depends so on one parameter in the equation, the dimensionless Λ−1 defined in (2.1). There is one special value Λ−1 c −1 is near Λ−1 that we have a codimension 2 bifurcation for Λ−1 = Λ−1 c , there is competition between c ; if Λ multiple modes, which leads to the following: (1) for Λ−1 < Λ−1 c , the equation (1.3) has a simple steady state bifurcation as we increase the bifurcation parameter σ; (2) for Λ−1 > Λ−1 c we will encounter a Hopf bifurcation followed by secondary bifurcations and finally we will reach a standing wave solution for σ sufficiently large. We also observe an interesting phenomenon: the amplitude of alternans at the stimulus site is quite different for the steady state and the periodic solution. For instance, the pure z-mode (standing wave solution, as shown in Fig. 3d) has an nonzero constant value at the stimulus site, but the pure r-mode solution (pure periodic, as shown in Fig. 3b) almost has no amplitude of oscillation at x = 0. This can be seen from the eigenfunctions φ′j s (see Fig. 5A). The z-mode is approximately given by φ0 and the r-mode, by φ1,2 . The different behavior between φ0 and φ1,2 induces the great change of the amplitude of solution at x = 0. Throughout this paper, the dimensionless cardiac fiber length is assumed to be ℓ̄ = 6. However, the dynamical behavior of (1.3) is similar for any ℓ̄ greater than 6. In particular, the evolution of eigenvalues has the same behavior as in Figure 1 when we increase Λ−1 . For smaller ℓ̄ the behavior of (1.3) is different, as we shall explore in a future paper. 1238 S. DAI AND D.G. SCHAEFFER σ = 12.0 σ = 12.4 σ = 12.8 0.8 0.8 0.8 0.4 0.4 0.4 0 0 0 −0.4 −0.4 −0.4 −0.8 −0.8 −0.4 0 A 0.4 0.8 −0.8 −0.8 −0.4 0 0.4 B 0.8 −0.8 −0.8 −0.4 0 0.4 0.8 C Figure 8. Phase plane portraits of the solutions for various values of σ, assuming ℓ̄ = 15 and Λ̄−1 = 10, which is obtained by plotting a(y2 , t) versus a(y1 , t), where y1 = 0.75ℓ̄ and y2 = 0.85ℓ̄. Incidentally, for the Noble model, dimensionless length ℓ̄ = 6 corresponds to a physical length of 4.32 cm. This conclusion follows from undoing the scalings ℓ̄ = wℓ/ξ 2 and using w = 0.045 cm and ξ = 0.18 cm (see [3] for computation of w and ξ). −1 For Λ−1 in a broader neighborhood of Λ−1 c , possibly far away from Λc , (1.3) may have more complicated dynamics, including chaos, which we will investigate in a future paper. Figure 8 shows the phase plane portraits for three different values of σ with fixed ℓ̄ = 15 and Λ̄−1 = 10. It appears case A is periodic, B is period doubled and C is chaotic. Acknowledgements. We gratefully acknowledge the support of the National Institutes of Health under grant 1R01-HL72831 and the National Science Foundation under grant PHY-0549259. References [1] J. Carr, Applications of Centre Manifold Theory. Springer-Verlag, New York (1981). [2] S. Dai and D.G. Schaeffer, Spectrum of a linearized amplitude equation for alternans in a cardiac fiber. SIAM J. Appl. Math. 69 (2008) 704–719. [3] B. Echebarria and A. Karma, Instability and spatiotemporal dynamics of alternans in paced cardiac tissue. Phys. Rev. Lett. 88 (2002) 208101. [4] B. Echebarria and A. Karma, Amplitude-equation approach to spatiotemporal dynamics of cardiac alternans. Phys. Rev. E 76 (2007) 051911. [5] A. Garfinkel, Y.-H. Kim, O. Voroshilovsky, Z. Qu, J.R. Kil, M.-H. Lee, H.S. Karagueuzian, J.N. Weiss and P.-S. Chen, Preventing ventricular fibrillation by flattening cardiac restitution. Proc. Natl. Acad. Sci. USA 97 (2000) 6061–6066. [6] R.F. Gilmour Jr. and D.R. Chialvo, Electrical restitution, Critical mass, and the riddle of fibrillation. J. Cardiovasc. Electrophysiol. 10 (1999) 1087–1089. [7] M. Golubitsky and D.G. Schaeffer, Singularities and Groups in Bifurcation Theory. Springer-Verlag, New York (1985). [8] J. Guckenheimer, On a codimension two bifurcation, in Dynamical Systems and Turbulence, Warwick 1980, Lect. Notes in Mathematics 898, Springer (1981) 99–142. [9] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dyanamical Systems, and Bifurcations of Vector Fields. SpringerVerlag, New York (1983). [10] M.R. Guevara, G. Ward, A. Shrier and L. Glass, Electrical alternans and period doubling bifurcations, in Proceedings of the 11th Computers in Cardiology Conference, IEEE Computer Society, Los Angeles, USA (1984) 167–170. [11] P. Holmes, Unfolding a degenerate nonlinear oscillator: a codimension two bifurcation, in Nonlinear Dynamics, R.H.G. Helleman Ed., New York Academy of Sciences, New York (1980) 473–488. [12] W.F. Langford, Periodic and steady state interactions lead to tori. SIAM J. Appl. Math. 37 (1979) 22–48. [13] C.C. Mitchell and D.G. Schaeffer, A two-current model for the dynamics of the cardiac membrane. Bull. Math. Biol. 65 (2003) 767–793. [14] D. Noble, A modification of the Hodgkin-Huxley equations applicable to Purkinje fiber actoin and pacemaker potential. J. Physiol. 160 (1962) 317–352. [15] J.B. Nolasco and R.W. Dahlen, A graphic method for the study of alternation in cardiac action potentials. J. Appl. Physiol. 25 (1968) 191–196. [16] A.V. Panfilov, Spiral breakup as a model of ventricular fibrillation. Chaos 8 (1998) 57–64.