- Email: [email protected]
A kinematic study of spiral wave drift due to an electric field
Physics Letters A 308 (2003) 179–186 www.elsevier.com/locate/pla
A kinematic study of spiral wave drift due to an electric field Zengru Di a,b , Zhilin Qu a,∗ , James N. Weiss a , Alan Garfinkel a a Department of Medicine (Cardiology), University of California, Los Angeles, CA 90095, USA b Department of Physics, Beijing Normal University, Beijing 100875, PR China
Received 8 July 2002; received in revised form 26 December 2002; accepted 2 January 2003 Communicated by A.R. Bishop
Abstract Kinematic equations are developed to describe the drift of a spiral wave under an external electric field. We find that the kinematic equations explain several observed properties of spiral wave drift. In particular, it explains the observation of Krinsky et al. Phys. Rev. Lett. 76 (1996) 3854 that dense spiral waves drift oppositely to sparse spiral waves. 2003 Elsevier Science B.V. All rights reserved. PACS: 82.40.Ck; 87.19. Hh; 05.45.-a; 47.32.Cc
Spiral wave drift induced by electric fields in excitable media have been widely studied [1–5]. It is well known that drift can occur not only along the electric field but also perpendicular to it. In a recent study, Krinsky et al. [1] showed that dense and sparse spiral waves drifted in opposite directions under the same electric field. In this Letter, we provide a kinematical theory [6–8] of spiral wave drift, which accounts for these observations. We consider an excitable medium with an electric field applied parallel to the x-axis and in which only one variable diffuses, with a unit diffusion constant [1,2], ∂u/∂t = f (u, v) + E∂u/∂x + ∇ 2 u, ∂v/∂t = g(u, v) + δE∂v/∂x, * Corresponding author.
E-mail address: [email protected] (Z. Qu).
(1)
where E is the strength of the electric field and δ is a coefficient. By changing to a moving coordinate frame x = x + δEt, the v convection term in the second equation of Eq. (1) can always be removed, as if δ = 0. If δ = 1 both convection terms in Eq. (1) can be removed, indicating that a spiral wave in Eq. (1) with δ = 1 must drift with a velocity −E along the direction of the electric field. Therefore, we only need to deal with the case of δ = 1, and the general case, for which we can consider δ = 0. In Eq. (1), we use Barkley’s kinetics [9]: f (u, v) = u(1 − u)[u − (v + b)/a]/ε, and g(u, v) = u − v. By varying a, b or ε, one can generate either dense or sparse spiral waves which drift in opposite directions [1]. In numerical simulation of Eq. (1), we used a five-point Euler method with a time step of 0.0006 time unit and a space step of 0.15625 (= 40/256) space unit. The domain size is 40 × 40. The same time step and space step were used in simulating Eq. (2) below. There are several morphologically significant features of a stable spiral wave (Fig. 1A, also see
0375-9601/03/$ – see front matter 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0375-9601(03)00031-8
180
Z. Di et al. / Physics Letters A 308 (2003) 179–186
Fig. 1. A schematic diagram of a stable spiral wave. O, Q, and q are detailed in the text. cn —normal velocity of Q; cg —tangential velocity of Q. θ —angle between OQ and the x-axis. Dashed circle is the trace of Q. B κ versus θ for the spiral wave. C cn versus κ from the 2D spiral wave. The parameters for the Barkley model are a = 1.0, b = 0.05, and ε = 0.02, which are in the dense spiral wave regime.
Ref. [6]): the rotation center O; the point q where wavefront meets waveback; and the point Q where the straight line from O to Q is tangent to the wavefront. The normal velocity (cn ) of Q and the spiral wave period Ts are related by cn = ωRQ = 2πRQ /Ts , where ω is the angular velocity and RQ the distance from O to Q. Following Zykov [6], Q is the point that has a critical curvature, beyond which conduction fails. In a spiral wave, curvature is low in the distal arm but increases inwardly to a critical curvature at Q. In other words, Q is the point distinguishing active (or regenerative) propagation and passive (or decremental) propagation [6]. We numerically solved Eq. (1) in a stable spiral wave parameter regime and determined O, Q, and q geometrically. In determining Q and q and calculating the curvature of the spiral wave, we used the contour line of u = 0.4 as the spiral wave shown in Fig. 1A. Fig. 1B shows the calculated the curvature (κ) of the wavefront and the waveback for a dense spiral wave, showing that curvature increases from the distal arm to a maximum (∼ 1.2) at Q, then decreases from Q to zero at q. We also calculated cn for the wavefront for this same spiral wave; it decreases from about 3 to 1.2 at Q (solid line in Fig. 1C), and then to zero at q (dashed line in Fig. 1C). Note that although from Q to q, curvature decreases from the maximum to zero, the normal velocity also decreases to zero. To study the effects of curvature and electric field on propagation, we numerically simulated the following equations {see Eqs. (6.4), (6.5) in Ref. [6] and Ref. [10]} in a 1-dimensional ring: ∂u/∂t = f (u, v) + Ez ∂u/∂z + κ∂u/∂z + ∂ 2 u/∂z2 ,
∂v/∂t = g(u, v) + δEz ∂v/∂z,
(2)
where Ez is the electric field strength along the ring. We initiated a unidirectional reentry wave in a long ring and shortened it gradually until conduction failed. For each shortening, we let the reentrant wave propagate for at least 10 cycles to reach its stationary state, and then we measured the cycle time T and conduction velocity c. Fig. 2A shows a family of such curves with parameters the same as for the dense spiral wave in Fig. 1 and Ez = 0. The last point of each curve is the point at which conduction fails. These points form a line, which we denote as cf = η(κ) or Tf = σ (κ). To verify that Q is the point which has the critical curvature at which conduction fails, we draw a vertical line (dashed line in Fig. 2A) at the spiral wave period Ts = 3.5. It intersects with the conduction failure line at Q . At Q , we have c = 1.2 and κ = 1.1, which are virtually identical to the values at the Q point measured in the spiral wave (Fig. 1C). We also did the same calculations for the sparse spiral wave and show the dispersion curves in Fig. 2B. In this case, the spiral wave period is Ts = 5.25 and the geometrically determined cn for the Q point is 2.1 (we also used the contour line of u = 0.4 as for the dense spiral wave). The vertical line of Ts = 5.25 intersects the conduction failure line at Q with c = 2 and κ = 1.06, which agrees with the geometrically determined cn for the Q point. The results shown in Fig. 2 indicate that the Q point is the critical point distinguishing active and decremental conduction. Its normal velocity cn is the cf in the ring for a κ similar to the curvature at the Q point of the spiral wave. We observed a similar feature in an earlier study in a cardiac tissue
Z. Di et al. / Physics Letters A 308 (2003) 179–186
181
Fig. 2. Dispersion curves for different curvatures. A dense spiral wave regime with a = 1.0, b = 0.05, and ε = 0.02. B sparse spiral wave regime with a = 0.55, b = 0.03, and ε = 0.02. Conduction failure is marked by the solid circles and the thick solid line. The dashed vertical line marks the spiral wave period Ts , which intersects the conduction failure line at Q .
model [11]. Note that the variability of c is larger for the case of dense spiral wave than the case of sparse spiral wave. For example, for κ = 0, c changes from 4.5 to 1.6 as T decrease until conduction fails for the case of dense spiral wave (Fig. 2A), while c changes from 3.65 to 2.5 for the case of sparse spiral wave (Fig. 2B). As κ increases, this variability is markedly reduced in the case of sparse spiral wave. It should be pointed out that the curvature for the point Q in the dense spiral wave happens to be at its maximum, but this is not the case in the sparse spiral wave. When Ez = 0 in Eq. (2), cf and Tf can be expressed as: cf = η(κ) − Ez ,
Tf = 2πRf /cf = 2πRf / η(κ) − Ez ,
for δ = 1, (3)
and cf = η(κ + Ez ), Tf = σ (κ + Ez ),
for δ = 0.
(4)
We calculated cf and Tf by simulating Eq. (2) for δ = 1 and δ = 0 (Fig. 2). We calculated cf and Tf for three types of drift: SW1—a dense spiral wave that drifts normally (Fig. 3A and D); SW2—a spiral wave that drifts perpendicular to the electrical field (Fig. 3B and E); SW3—a sparse spiral wave that drifts oppositely to SW1 (Fig. 3C and F).
We now turn to formulating our kinematic equations for the spiral wave drift under the electric field. Because Q is the critical point that is extremely important for the spiral wave movement, we formulate differential equations to describe the motion of Q. Since the electric field changes the movement of the spiral wave, a geometrical definition of Q as in Fig. 1A is not applicable. Here, we generally define the Q point as the point that distinguishes active and passive propagation in the spiral wave. At point Q, there is a normal velocity cn and a tangential velocity cg (Fig. 1A). We assume the rotation frequency is ω. Then the differential equations for the motion of Q are: x˙ = cg cos θ − cn sin θ, y˙ = cg sin θ + cn cos θ, θ˙ = ω,
(5)
where x and y are the coordinates of Q and θ is the angle between OQ and the x-axis as, shown in Fig. 1A. Once we know cn , cg and ω, we can solve Eq. (5) to get the trajectory of Q. For a stable spiral wave in Eq. (1) without an electric field, the velocity of Q has only its normal component, i.e., cn = η(κ0 ), cg = 0, and a fixed angular velocity, ω = 2π/Ts . κ0 is the critical curvature at Q, and Ts is the period of the spiral wave. Eq. (5) then result in a circle for Q.
182
Z. Di et al. / Physics Letters A 308 (2003) 179–186
Fig. 3. cf and Tf versus Ez for a = 1.0, b = 0.05, and ε = 0.02 A–D; a = 0.75, b = 0.05, ε = 0.025 B–E; and a = 0.55, b = 0.03, and ε = 0.02 C–F for δ = 1 (solid circles and lines) and δ = 0 (open circles and dashed lines). We used, in Eq. (2), κ = 1.1 for A–D and 1.06 for E, F. These curvatures were taken from Q as in Fig. 2.
For a spiral wave under an electric field with δ = 1 in Eq. (1), the velocity of Q changes. We decompose the electric field E at Q into a normal component E sin θ and a tangential component −E cos θ . Following Eq. (5), we assume that the normal velocity is cn = η(κ0 ) + E sin θ , and, following that, cg = −E cos θ . Considering that in the moving frame of velocity −E, the spiral wave has period Ts and cn = η(κ0 ), we can set ω = 2π/Ts . With this setting of cn , cg , and ω, Eq. (5) results in a spiral wave that drifts along the direction of the electric field with a velocity −E. For a spiral wave under an electric field with δ = 0, we choose cn as in Eq. (4), i.e.,
cn = η(κ0 − E sin θ ),
(6)
but cg and ω cannot be explicitly determined as for δ = 1. We phenomenologically formulate cg and ω according to some critical observations. One observation in both simulations and experiments is that for δ = 1, there is no deformation of the spiral wave, but for δ = 1 , deformation of spiral wave occurs [2,3,5]. This deformation results in a curvature change in the spiral wave, which then causes a tangential movement of the Q point, because the spiral wave has to adjust Q to the point that has the critical curvature. Another observation is that for δ = 1, cn changes with E sin θ
Z. Di et al. / Physics Letters A 308 (2003) 179–186
nonlinearly, and so does Tf , which forces Q to adjust its position constantly when θ changes. This will also induce a tangential velocity to Q. Therefore, we correct cg from the case of δ = 1 as follows: cg = −E cos θ + "cg ,
(7)
where "cg = α η(κ0 − E sin θ ) − η(κ0 ) + E sin θ + ρ σ (κ0 − E sin θ ) × η(κ0 − E sin θ ) − 2πRQ .
(8)
The first term in the right-hand side of Eq. (8) is due to the deformation of the spiral wave, which we assume is proportional to the normal velocity difference between the δ = 0 case and the δ = 1 case. The second term is due to the nonlinear change of cn and Tf caused by the electrical field, which we assume that it is proportional to the difference of the product of cn and Tf between the δ = 0 case and the δ = 1 case. α and ρ in Eq. (8) are two adjustable parameters. We assume ω changes as follows: ω = 2π/(Ts + γ "T ) − (β"cg /cn ), (9) where "T = σ (κ − E sin θ ) − 2πRQ / η(κ0 ) + E sin θ , (10) which is the difference in Tf between the δ = 0 case and the δ = 1 case. The second correction term in Eq. (9) is caused by the correction "cg to cg , and β and γ are adjustable parameters. Having cn , cg and ω, we can simulate Eq. (5) and compare it to results obtained from a 2D spiral wave in the electric field. Fig. 4A–C show the traces of Q from the simulation of Eq. (5) and tip trajectories from the 2D simulation of Eq. (1). Functions η and σ were fitted from data shown Fig. 2.1 In Fig. 4D–F, 1 Functions η and σ fitted from Fig. 2 for the three types of spiral wave. SW1: η(κ) = 1.56 − 0.424κ + 0.0822κ 2 , σ (κ) = 0.305 + 3.21 exp[(κ − 1.21)/2.28] + 0.0001 exp[(κ − 1.21)/0.178], and κ0 = 1.1. SW2: η(κ) = 1.96 − 0.473κ + 0.0587κ 2 , σ (κ) = 0.763 + 2.78 exp[(κ − 0.772)/1.94] + 0.00776 exp[(κ − 0.772)/0.156], and κ0 = 1.1. SW3: η(κ) = 1.56 − 0.424κ + 0.0822κ 2 , σ (κ) = 0.887 + 2.39 exp[(κ − 0.446)/2.17] + 0.0696 exp[(κ − 0.446)/0.218], if κ > 1.15, then η(κ) = η(1.15), and σ (κ) = σ (1.15), κ0 = 1.06. We used second order polynomial to fit the velocity and double exponential
183
we show the drifting velocities in the x-direction (vx ) and in the y-direction (vy ) for the three cases, respectively. Note that the only change in each case is the electric field E (from 0 to 0.3), no other parameter changes in either the 2D simulation of Eq. (1) or the kinematic simulation of Eq. (5). This indicates that the functions η and σ are crucial for determining the drifting properties. In other words, the drift of the spiral wave is highly determined by the propagation properties shown Fig. 3, though we have four free parameters to choose. To understand how the rotation changes, we plot the tip trajectory in the moving frame so that it forms a closed loop. We then calculated the tip velocity in this frame. Fig. 5 shows the tip trajectories (left) and the tip velocities (right) from both the 2D spiral wave simulation of Eq. (1) and the kinematic simulation of Eq. (5). For δ = 1, the tip trajectory is a circle and the tip velocity is constant, as expected. For δ = 0, the tip trajectories are no longer a circle, but are somewhat deformed. The tip velocities are also not constant (Fig. 5D and F), indicating that the instantaneous rotation frequency is also changing. This supports our correction to ω. Krinsky et al. [1] suggested two mechanisms for spiral wave drift under an electric field: periodic modulation of the core size and periodic modulation of velocity. The competition between these two mechanisms results in different drifting directions for dense and sparse spiral waves. Our kinematic study supports their conclusion, but we have a more specific explanation. Using our kinematic simulation, we have identified the following effects: (1) If we only correct cn as in Eq. (6), but not cg and ω, the spiral wave drifts slower than without this correction. But this effect is not strong enough to cause the spiral wave to drift in the opposite direction; (2) If we only correct cg as in Eqs. (7) and (8), but keep cn as η(κ0 ) + E sin θ and ω = 2π/Ts , it results in a perpendicular drift while the drift velocity remains −E;
for the period. Other choices of fitting function do not change our results. In Fig. 2, the data was plotted as functions of Ez , but according to Eq. (4), when δ = 0 κ and Ez are equivalent.
184
Z. Di et al. / Physics Letters A 308 (2003) 179–186
Fig. 4. A–C tip trajectories for SW1(A), SW2(B), and SW3(C). Thick lines are from a 2D spiral wave simulation, thin lines are from our kinematic simulation. The solid arrow indicates the direction of the electric field. Solid arrows indicate the drift directions. In 2D spiral wave simulation, we used the intersection of two consecutive contour lines of u = 0.4 as the tip. This intersection lies between Q and q, which resulted in a smaller trajectory from the 2D simulation than from the kinematic simulation. This is also the cause that the velocity shown in Fig. 5 is smaller for 2D simulation than for kinematic simulation. D–F drift velocity vx (circles) and vy (squares) versus E for SW1(D), SW2(E), and SW3(F). Solid symbols are from the 2D spiral wave simulation, and open symbols are from the kinematic simulation. We used α = 1.24, β = 0.7, γ = 0.85, and ρ = 1.2 for SW1; α = 3.5, β = 0.65, γ = 1.2, and ρ = 3.0 for SW2; and α = 1.2, β = 1.0, γ = 1.3, and ρ = 1.2 for SW3. E = 0.2 and δ = 0. In simulation, we first fitted the tip trajectory for a certain E (as in A–C) by adjusting α, β, γ , and ρ, and then change E to see whether this parameter set can fit the drifting velocities for all E.
(3) If we only correct ω, as in Eq. (9), it can induce a drift either parallel or anti-parallel to the electric field. For a dense spiral wave, it induces an antiparallel drift, while for a sparse spiral wave induces a parallel drift. The reason is that the corrections in "T for dense and sparse spiral wave have different signs, according to Eq. (10) and Fig. 3. This is the major qualitative difference between the dense and the sparse spiral wave.
Therefore, perpendicular drift is induced by nonlinear corrections in cg , and parallel and antiparallel drift are induced by corrections in cn and ω. For dense and sparse spiral waves, the corrections in ω are of opposite sign, which causes the opposite drift. In conclusion, we formulated a set of kinematic equations to describe the spiral wave drift under an external electric field, and incorporated the possible factors that may regulate the tip motion of the spiral
Z. Di et al. / Physics Letters A 308 (2003) 179–186
185
Fig. 5. Tip trajectories (left panels) and velocities (right panels) in the moving frames. In the left panels, thick lines are from 2D spiral wave simulations and thin lines are from our kinematic simulations of Eq. (5). In the right panels, tip velocity (thick lines) and their x (solid thin lines) and y (dashed lines) components are from a 2D spiral wave (upper panel) and from the kinematic simulation (lower panel). E = 0.2. A SW3 for δ = 1. B SW1 for δ = 0. C SW3 for δ = 0.
wave. We used conduction velocity and rotation period at conduction failure from the simulation of a 1D ring as functions of curvature or strength of the electric field. For each spiral wave, our kinematic equations simulated qualitatively and quantitatively the drift of the spiral wave for different strength of the electric field. Our study shows that the conduction properties at conduction failure of reentry determine the properties of spiral wave under electric field. In a recent study [4], Hakim and Karma showed that spiral waves in weakly excitable media drift along the direction of an electric field, opposite to the normal
case. We anticipate that spiral waves in ischemic cardiac tissue (where excitability is low) will drift oppositely to spiral waves in normal tissue.
Acknowledgements This research was supported by NIH SCOR in Sudden Cardiac Death P50 HL52319, by American Heart Association (Z.Q.), and by Education Council of People’s Republic of China (Z.D.).
186
Z. Di et al. / Physics Letters A 308 (2003) 179–186
References [1] V.I. Krinsky, E. Hamm, V. Voignier, Phys. Rev. Lett. 76 (1996) 3854. [2] A.P. Munuzuri, M. Gomez-Gesteira, V. Perez-Munuzuri, et al., Phys. Rev. E 48 (1993) R3232. [3] O. Steinbock, J. Schutze, C. Muller, Phys. Rev. Lett. 68 (1992) 248. [4] V. Hakim, A. Karma, Phys. Rev. E 60 (1999) 5073.
[5] M. Gabbay, E. Ott, P.N. Guzdar, Phys. Rev. E 59 (1999) 2443. [6] V.S. Zykov, Simulation of Wave Processes in Excitable Media, Manchester Univ. Press, Manchester, 1982. [7] A.S. Mikhailov, V.S. Zykov, Physica D 52 (1991) 379. [8] A.S. Mikhailov, V.A. Davydov, V.S. Zykov, Physica D 70 (1994) 1. [9] D. Barkley, Phys. Rev. Lett. 68 (1992) 2090. [10] P. Comtois, A. Vinet, Phys. Rev. E 60 (1999) 4619. [11] Z. Qu, F. Xie, A. Garfinkel, et al., Ann. Biomed. Engrg. 28 (2000) 755.
A kinematic study of spiral wave drift due to an electric field Zengru Di a,b , Zhilin Qu a,∗ , James N. Weiss a , Alan Garfinkel a a Department of Medicine (Cardiology), University of California, Los Angeles, CA 90095, USA b Department of Physics, Beijing Normal University, Beijing 100875, PR China
Received 8 July 2002; received in revised form 26 December 2002; accepted 2 January 2003 Communicated by A.R. Bishop
Abstract Kinematic equations are developed to describe the drift of a spiral wave under an external electric field. We find that the kinematic equations explain several observed properties of spiral wave drift. In particular, it explains the observation of Krinsky et al. Phys. Rev. Lett. 76 (1996) 3854 that dense spiral waves drift oppositely to sparse spiral waves. 2003 Elsevier Science B.V. All rights reserved. PACS: 82.40.Ck; 87.19. Hh; 05.45.-a; 47.32.Cc
Spiral wave drift induced by electric fields in excitable media have been widely studied [1–5]. It is well known that drift can occur not only along the electric field but also perpendicular to it. In a recent study, Krinsky et al. [1] showed that dense and sparse spiral waves drifted in opposite directions under the same electric field. In this Letter, we provide a kinematical theory [6–8] of spiral wave drift, which accounts for these observations. We consider an excitable medium with an electric field applied parallel to the x-axis and in which only one variable diffuses, with a unit diffusion constant [1,2], ∂u/∂t = f (u, v) + E∂u/∂x + ∇ 2 u, ∂v/∂t = g(u, v) + δE∂v/∂x, * Corresponding author.
E-mail address: [email protected] (Z. Qu).
(1)
where E is the strength of the electric field and δ is a coefficient. By changing to a moving coordinate frame x = x + δEt, the v convection term in the second equation of Eq. (1) can always be removed, as if δ = 0. If δ = 1 both convection terms in Eq. (1) can be removed, indicating that a spiral wave in Eq. (1) with δ = 1 must drift with a velocity −E along the direction of the electric field. Therefore, we only need to deal with the case of δ = 1, and the general case, for which we can consider δ = 0. In Eq. (1), we use Barkley’s kinetics [9]: f (u, v) = u(1 − u)[u − (v + b)/a]/ε, and g(u, v) = u − v. By varying a, b or ε, one can generate either dense or sparse spiral waves which drift in opposite directions [1]. In numerical simulation of Eq. (1), we used a five-point Euler method with a time step of 0.0006 time unit and a space step of 0.15625 (= 40/256) space unit. The domain size is 40 × 40. The same time step and space step were used in simulating Eq. (2) below. There are several morphologically significant features of a stable spiral wave (Fig. 1A, also see
0375-9601/03/$ – see front matter 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0375-9601(03)00031-8
180
Z. Di et al. / Physics Letters A 308 (2003) 179–186
Fig. 1. A schematic diagram of a stable spiral wave. O, Q, and q are detailed in the text. cn —normal velocity of Q; cg —tangential velocity of Q. θ —angle between OQ and the x-axis. Dashed circle is the trace of Q. B κ versus θ for the spiral wave. C cn versus κ from the 2D spiral wave. The parameters for the Barkley model are a = 1.0, b = 0.05, and ε = 0.02, which are in the dense spiral wave regime.
Ref. [6]): the rotation center O; the point q where wavefront meets waveback; and the point Q where the straight line from O to Q is tangent to the wavefront. The normal velocity (cn ) of Q and the spiral wave period Ts are related by cn = ωRQ = 2πRQ /Ts , where ω is the angular velocity and RQ the distance from O to Q. Following Zykov [6], Q is the point that has a critical curvature, beyond which conduction fails. In a spiral wave, curvature is low in the distal arm but increases inwardly to a critical curvature at Q. In other words, Q is the point distinguishing active (or regenerative) propagation and passive (or decremental) propagation [6]. We numerically solved Eq. (1) in a stable spiral wave parameter regime and determined O, Q, and q geometrically. In determining Q and q and calculating the curvature of the spiral wave, we used the contour line of u = 0.4 as the spiral wave shown in Fig. 1A. Fig. 1B shows the calculated the curvature (κ) of the wavefront and the waveback for a dense spiral wave, showing that curvature increases from the distal arm to a maximum (∼ 1.2) at Q, then decreases from Q to zero at q. We also calculated cn for the wavefront for this same spiral wave; it decreases from about 3 to 1.2 at Q (solid line in Fig. 1C), and then to zero at q (dashed line in Fig. 1C). Note that although from Q to q, curvature decreases from the maximum to zero, the normal velocity also decreases to zero. To study the effects of curvature and electric field on propagation, we numerically simulated the following equations {see Eqs. (6.4), (6.5) in Ref. [6] and Ref. [10]} in a 1-dimensional ring: ∂u/∂t = f (u, v) + Ez ∂u/∂z + κ∂u/∂z + ∂ 2 u/∂z2 ,
∂v/∂t = g(u, v) + δEz ∂v/∂z,
(2)
where Ez is the electric field strength along the ring. We initiated a unidirectional reentry wave in a long ring and shortened it gradually until conduction failed. For each shortening, we let the reentrant wave propagate for at least 10 cycles to reach its stationary state, and then we measured the cycle time T and conduction velocity c. Fig. 2A shows a family of such curves with parameters the same as for the dense spiral wave in Fig. 1 and Ez = 0. The last point of each curve is the point at which conduction fails. These points form a line, which we denote as cf = η(κ) or Tf = σ (κ). To verify that Q is the point which has the critical curvature at which conduction fails, we draw a vertical line (dashed line in Fig. 2A) at the spiral wave period Ts = 3.5. It intersects with the conduction failure line at Q . At Q , we have c = 1.2 and κ = 1.1, which are virtually identical to the values at the Q point measured in the spiral wave (Fig. 1C). We also did the same calculations for the sparse spiral wave and show the dispersion curves in Fig. 2B. In this case, the spiral wave period is Ts = 5.25 and the geometrically determined cn for the Q point is 2.1 (we also used the contour line of u = 0.4 as for the dense spiral wave). The vertical line of Ts = 5.25 intersects the conduction failure line at Q with c = 2 and κ = 1.06, which agrees with the geometrically determined cn for the Q point. The results shown in Fig. 2 indicate that the Q point is the critical point distinguishing active and decremental conduction. Its normal velocity cn is the cf in the ring for a κ similar to the curvature at the Q point of the spiral wave. We observed a similar feature in an earlier study in a cardiac tissue
Z. Di et al. / Physics Letters A 308 (2003) 179–186
181
Fig. 2. Dispersion curves for different curvatures. A dense spiral wave regime with a = 1.0, b = 0.05, and ε = 0.02. B sparse spiral wave regime with a = 0.55, b = 0.03, and ε = 0.02. Conduction failure is marked by the solid circles and the thick solid line. The dashed vertical line marks the spiral wave period Ts , which intersects the conduction failure line at Q .
model [11]. Note that the variability of c is larger for the case of dense spiral wave than the case of sparse spiral wave. For example, for κ = 0, c changes from 4.5 to 1.6 as T decrease until conduction fails for the case of dense spiral wave (Fig. 2A), while c changes from 3.65 to 2.5 for the case of sparse spiral wave (Fig. 2B). As κ increases, this variability is markedly reduced in the case of sparse spiral wave. It should be pointed out that the curvature for the point Q in the dense spiral wave happens to be at its maximum, but this is not the case in the sparse spiral wave. When Ez = 0 in Eq. (2), cf and Tf can be expressed as: cf = η(κ) − Ez ,
Tf = 2πRf /cf = 2πRf / η(κ) − Ez ,
for δ = 1, (3)
and cf = η(κ + Ez ), Tf = σ (κ + Ez ),
for δ = 0.
(4)
We calculated cf and Tf by simulating Eq. (2) for δ = 1 and δ = 0 (Fig. 2). We calculated cf and Tf for three types of drift: SW1—a dense spiral wave that drifts normally (Fig. 3A and D); SW2—a spiral wave that drifts perpendicular to the electrical field (Fig. 3B and E); SW3—a sparse spiral wave that drifts oppositely to SW1 (Fig. 3C and F).
We now turn to formulating our kinematic equations for the spiral wave drift under the electric field. Because Q is the critical point that is extremely important for the spiral wave movement, we formulate differential equations to describe the motion of Q. Since the electric field changes the movement of the spiral wave, a geometrical definition of Q as in Fig. 1A is not applicable. Here, we generally define the Q point as the point that distinguishes active and passive propagation in the spiral wave. At point Q, there is a normal velocity cn and a tangential velocity cg (Fig. 1A). We assume the rotation frequency is ω. Then the differential equations for the motion of Q are: x˙ = cg cos θ − cn sin θ, y˙ = cg sin θ + cn cos θ, θ˙ = ω,
(5)
where x and y are the coordinates of Q and θ is the angle between OQ and the x-axis as, shown in Fig. 1A. Once we know cn , cg and ω, we can solve Eq. (5) to get the trajectory of Q. For a stable spiral wave in Eq. (1) without an electric field, the velocity of Q has only its normal component, i.e., cn = η(κ0 ), cg = 0, and a fixed angular velocity, ω = 2π/Ts . κ0 is the critical curvature at Q, and Ts is the period of the spiral wave. Eq. (5) then result in a circle for Q.
182
Z. Di et al. / Physics Letters A 308 (2003) 179–186
Fig. 3. cf and Tf versus Ez for a = 1.0, b = 0.05, and ε = 0.02 A–D; a = 0.75, b = 0.05, ε = 0.025 B–E; and a = 0.55, b = 0.03, and ε = 0.02 C–F for δ = 1 (solid circles and lines) and δ = 0 (open circles and dashed lines). We used, in Eq. (2), κ = 1.1 for A–D and 1.06 for E, F. These curvatures were taken from Q as in Fig. 2.
For a spiral wave under an electric field with δ = 1 in Eq. (1), the velocity of Q changes. We decompose the electric field E at Q into a normal component E sin θ and a tangential component −E cos θ . Following Eq. (5), we assume that the normal velocity is cn = η(κ0 ) + E sin θ , and, following that, cg = −E cos θ . Considering that in the moving frame of velocity −E, the spiral wave has period Ts and cn = η(κ0 ), we can set ω = 2π/Ts . With this setting of cn , cg , and ω, Eq. (5) results in a spiral wave that drifts along the direction of the electric field with a velocity −E. For a spiral wave under an electric field with δ = 0, we choose cn as in Eq. (4), i.e.,
cn = η(κ0 − E sin θ ),
(6)
but cg and ω cannot be explicitly determined as for δ = 1. We phenomenologically formulate cg and ω according to some critical observations. One observation in both simulations and experiments is that for δ = 1, there is no deformation of the spiral wave, but for δ = 1 , deformation of spiral wave occurs [2,3,5]. This deformation results in a curvature change in the spiral wave, which then causes a tangential movement of the Q point, because the spiral wave has to adjust Q to the point that has the critical curvature. Another observation is that for δ = 1, cn changes with E sin θ
Z. Di et al. / Physics Letters A 308 (2003) 179–186
nonlinearly, and so does Tf , which forces Q to adjust its position constantly when θ changes. This will also induce a tangential velocity to Q. Therefore, we correct cg from the case of δ = 1 as follows: cg = −E cos θ + "cg ,
(7)
where "cg = α η(κ0 − E sin θ ) − η(κ0 ) + E sin θ + ρ σ (κ0 − E sin θ ) × η(κ0 − E sin θ ) − 2πRQ .
(8)
The first term in the right-hand side of Eq. (8) is due to the deformation of the spiral wave, which we assume is proportional to the normal velocity difference between the δ = 0 case and the δ = 1 case. The second term is due to the nonlinear change of cn and Tf caused by the electrical field, which we assume that it is proportional to the difference of the product of cn and Tf between the δ = 0 case and the δ = 1 case. α and ρ in Eq. (8) are two adjustable parameters. We assume ω changes as follows: ω = 2π/(Ts + γ "T ) − (β"cg /cn ), (9) where "T = σ (κ − E sin θ ) − 2πRQ / η(κ0 ) + E sin θ , (10) which is the difference in Tf between the δ = 0 case and the δ = 1 case. The second correction term in Eq. (9) is caused by the correction "cg to cg , and β and γ are adjustable parameters. Having cn , cg and ω, we can simulate Eq. (5) and compare it to results obtained from a 2D spiral wave in the electric field. Fig. 4A–C show the traces of Q from the simulation of Eq. (5) and tip trajectories from the 2D simulation of Eq. (1). Functions η and σ were fitted from data shown Fig. 2.1 In Fig. 4D–F, 1 Functions η and σ fitted from Fig. 2 for the three types of spiral wave. SW1: η(κ) = 1.56 − 0.424κ + 0.0822κ 2 , σ (κ) = 0.305 + 3.21 exp[(κ − 1.21)/2.28] + 0.0001 exp[(κ − 1.21)/0.178], and κ0 = 1.1. SW2: η(κ) = 1.96 − 0.473κ + 0.0587κ 2 , σ (κ) = 0.763 + 2.78 exp[(κ − 0.772)/1.94] + 0.00776 exp[(κ − 0.772)/0.156], and κ0 = 1.1. SW3: η(κ) = 1.56 − 0.424κ + 0.0822κ 2 , σ (κ) = 0.887 + 2.39 exp[(κ − 0.446)/2.17] + 0.0696 exp[(κ − 0.446)/0.218], if κ > 1.15, then η(κ) = η(1.15), and σ (κ) = σ (1.15), κ0 = 1.06. We used second order polynomial to fit the velocity and double exponential
183
we show the drifting velocities in the x-direction (vx ) and in the y-direction (vy ) for the three cases, respectively. Note that the only change in each case is the electric field E (from 0 to 0.3), no other parameter changes in either the 2D simulation of Eq. (1) or the kinematic simulation of Eq. (5). This indicates that the functions η and σ are crucial for determining the drifting properties. In other words, the drift of the spiral wave is highly determined by the propagation properties shown Fig. 3, though we have four free parameters to choose. To understand how the rotation changes, we plot the tip trajectory in the moving frame so that it forms a closed loop. We then calculated the tip velocity in this frame. Fig. 5 shows the tip trajectories (left) and the tip velocities (right) from both the 2D spiral wave simulation of Eq. (1) and the kinematic simulation of Eq. (5). For δ = 1, the tip trajectory is a circle and the tip velocity is constant, as expected. For δ = 0, the tip trajectories are no longer a circle, but are somewhat deformed. The tip velocities are also not constant (Fig. 5D and F), indicating that the instantaneous rotation frequency is also changing. This supports our correction to ω. Krinsky et al. [1] suggested two mechanisms for spiral wave drift under an electric field: periodic modulation of the core size and periodic modulation of velocity. The competition between these two mechanisms results in different drifting directions for dense and sparse spiral waves. Our kinematic study supports their conclusion, but we have a more specific explanation. Using our kinematic simulation, we have identified the following effects: (1) If we only correct cn as in Eq. (6), but not cg and ω, the spiral wave drifts slower than without this correction. But this effect is not strong enough to cause the spiral wave to drift in the opposite direction; (2) If we only correct cg as in Eqs. (7) and (8), but keep cn as η(κ0 ) + E sin θ and ω = 2π/Ts , it results in a perpendicular drift while the drift velocity remains −E;
for the period. Other choices of fitting function do not change our results. In Fig. 2, the data was plotted as functions of Ez , but according to Eq. (4), when δ = 0 κ and Ez are equivalent.
184
Z. Di et al. / Physics Letters A 308 (2003) 179–186
Fig. 4. A–C tip trajectories for SW1(A), SW2(B), and SW3(C). Thick lines are from a 2D spiral wave simulation, thin lines are from our kinematic simulation. The solid arrow indicates the direction of the electric field. Solid arrows indicate the drift directions. In 2D spiral wave simulation, we used the intersection of two consecutive contour lines of u = 0.4 as the tip. This intersection lies between Q and q, which resulted in a smaller trajectory from the 2D simulation than from the kinematic simulation. This is also the cause that the velocity shown in Fig. 5 is smaller for 2D simulation than for kinematic simulation. D–F drift velocity vx (circles) and vy (squares) versus E for SW1(D), SW2(E), and SW3(F). Solid symbols are from the 2D spiral wave simulation, and open symbols are from the kinematic simulation. We used α = 1.24, β = 0.7, γ = 0.85, and ρ = 1.2 for SW1; α = 3.5, β = 0.65, γ = 1.2, and ρ = 3.0 for SW2; and α = 1.2, β = 1.0, γ = 1.3, and ρ = 1.2 for SW3. E = 0.2 and δ = 0. In simulation, we first fitted the tip trajectory for a certain E (as in A–C) by adjusting α, β, γ , and ρ, and then change E to see whether this parameter set can fit the drifting velocities for all E.
(3) If we only correct ω, as in Eq. (9), it can induce a drift either parallel or anti-parallel to the electric field. For a dense spiral wave, it induces an antiparallel drift, while for a sparse spiral wave induces a parallel drift. The reason is that the corrections in "T for dense and sparse spiral wave have different signs, according to Eq. (10) and Fig. 3. This is the major qualitative difference between the dense and the sparse spiral wave.
Therefore, perpendicular drift is induced by nonlinear corrections in cg , and parallel and antiparallel drift are induced by corrections in cn and ω. For dense and sparse spiral waves, the corrections in ω are of opposite sign, which causes the opposite drift. In conclusion, we formulated a set of kinematic equations to describe the spiral wave drift under an external electric field, and incorporated the possible factors that may regulate the tip motion of the spiral
Z. Di et al. / Physics Letters A 308 (2003) 179–186
185
Fig. 5. Tip trajectories (left panels) and velocities (right panels) in the moving frames. In the left panels, thick lines are from 2D spiral wave simulations and thin lines are from our kinematic simulations of Eq. (5). In the right panels, tip velocity (thick lines) and their x (solid thin lines) and y (dashed lines) components are from a 2D spiral wave (upper panel) and from the kinematic simulation (lower panel). E = 0.2. A SW3 for δ = 1. B SW1 for δ = 0. C SW3 for δ = 0.
wave. We used conduction velocity and rotation period at conduction failure from the simulation of a 1D ring as functions of curvature or strength of the electric field. For each spiral wave, our kinematic equations simulated qualitatively and quantitatively the drift of the spiral wave for different strength of the electric field. Our study shows that the conduction properties at conduction failure of reentry determine the properties of spiral wave under electric field. In a recent study [4], Hakim and Karma showed that spiral waves in weakly excitable media drift along the direction of an electric field, opposite to the normal
case. We anticipate that spiral waves in ischemic cardiac tissue (where excitability is low) will drift oppositely to spiral waves in normal tissue.
Acknowledgements This research was supported by NIH SCOR in Sudden Cardiac Death P50 HL52319, by American Heart Association (Z.Q.), and by Education Council of People’s Republic of China (Z.D.).
186
Z. Di et al. / Physics Letters A 308 (2003) 179–186
References [1] V.I. Krinsky, E. Hamm, V. Voignier, Phys. Rev. Lett. 76 (1996) 3854. [2] A.P. Munuzuri, M. Gomez-Gesteira, V. Perez-Munuzuri, et al., Phys. Rev. E 48 (1993) R3232. [3] O. Steinbock, J. Schutze, C. Muller, Phys. Rev. Lett. 68 (1992) 248. [4] V. Hakim, A. Karma, Phys. Rev. E 60 (1999) 5073.
[5] M. Gabbay, E. Ott, P.N. Guzdar, Phys. Rev. E 59 (1999) 2443. [6] V.S. Zykov, Simulation of Wave Processes in Excitable Media, Manchester Univ. Press, Manchester, 1982. [7] A.S. Mikhailov, V.S. Zykov, Physica D 52 (1991) 379. [8] A.S. Mikhailov, V.A. Davydov, V.S. Zykov, Physica D 70 (1994) 1. [9] D. Barkley, Phys. Rev. Lett. 68 (1992) 2090. [10] P. Comtois, A. Vinet, Phys. Rev. E 60 (1999) 4619. [11] Z. Qu, F. Xie, A. Garfinkel, et al., Ann. Biomed. Engrg. 28 (2000) 755.