Spiral–pacemaker interactions in a mathematical model of excitable medium

Excitable media, such as flame fronts during a forest fire [1], the BZ reaction [2], and the heart [3] have certain features in common. They respond to a suprathreshold stimulus with a characteristic wave-like activity, e.g. a flame front, a chemical wave or an electrical wave. After the wave they cannot be excited again for some time, called the refractory period of the medium. The refractory period acts as the waveback of the excitation wave. Excitation waves propagate in space without any attenuation, but they annihilate each other upon collision. In two dimensional excitable media excitation waves are typically target waves emanating from a pacemaker or spiral waves. In spiral waves, the wavefront and waveback meet at the spiral tip. Winfree observed that the tip of spiral wave is not stationary but rotates around a circle (rigid rotation) or moves in complex flower like trajectories (meander or hyper meander) [4]. When the spiral meanders the activation near the spiral core is quasiperiodic [5–8]. Fluctuations in the spiral period because of this quasiperiodic activation is maximum near the core and diminishes as we move away from the core [9, 10]. Studies have shown that the transition from rigid rotation to meander as excitability of the system is varied is a supercritical Hopf bifurcation [11]. Meandering transition of a spiral wave is also related to the Euclidean symmetry of the spiral wave in a plane [11–13]. Apart from the rigid rotation and the meander, the spiral tip in the vicinity of heterogeneities can also drift from its initial location [3, 14].

Various types of nonexcitable or partially excitable heterogeneities in the medium can affect the dynamics of the spiral wave. For example experiments in cardiac tissue [15] and in cardiac cell cultures [16] have found that spiral waves interacting with a nonexcitable tissue or other obstacles can attach to them and assume a rotating frequency determined by the size of the obstacle. Detailed numerical work later established that such attachment depends on the local kinetics near the spiral tip [17] and size and position of the obstacle with respect to the spiral tip [18, 19]. Theoretical and experimental studies demonstrate that spiral waves trajectories can be strongly affected by partially excitable heterogeneities [20–23]. However, the dynamics can be more complex when, on occasion, the wave enters the heterogeneous region. Such complex dynamics were found in monolayers of cultured cardiac cells when the wave interacted with line defects and partial conduction blocks. Time series recorded from near the core of the spiral wave exhibited period-1, period-2 and period-3 oscillations [24, 25]. In numerical studies of wave dynamics in a medium with an ischemic region of reduced excitability Xu and Guevara [26] observed two kinds of spiral formation when the medium was paced with high frequency waves. In the first kind the spiral tip stayed entirely outside the ischemic zone and in the second kind the tip went around the ischemic zone but then returned to the original position through the heterogeneity.

Pumir et al [27] have designed an ischemic zone in monolayers of neonatal rat ventricular cells to study ectopic beat generation and spiral wave initiation. Similarly, we found a novel scenario in which spiral waves are initiated in cultured monolayers of chicken heart cells by the interaction of two pacemakers. We had engineered a pacemaker in the center of the monolayer and induced extra pacemakers in the periphery by applying potassium channel blockers to the monolayer [28]. The dynamic interaction between target waves from the faster side pacemaker in the medium with the central pacemaker led to wave breaks which subsequently formed spiral waves. The spiral wave thus generated then interacts with the pacemaker.

This observation motivates the current work in which we investigate spiral–pacemaker interactions. The spiral–pacemaker interactions can lead to complex spatiotemporal activation in the medium. The resulting dynamics depend on both on the period and the size of the pacemaker. We present the model equations in section 2, and our results in section 3. Implications of these results are discussed in section 4.

We use a modified version of FitzHugh–Nagumo (FHN) equations [28, 29], to simulate excitable media and pacemaker dynamics. The equations are:

$$\begin{eqnarray} \frac{\partial v}{\partial t} &=& \frac{1}{\varepsilon}(v-v^3/3-w)+I_{\mathrm{P}}+D\frac{\partial^2 v}{\partial x^2},\nonumber \\ \frac{\partial w}{\partial t} &=& (v+\beta-\gamma w)\left(\frac{w_{\mathrm{H}}-w_{\mathrm{L}}}{1+\mathrm{e}^{-4v}}+w_{\mathrm{L}}\right), \end{eqnarray} \tag{ 1 }$$

where v is the excitation variable and w is the recovery variable. The model equation is modified from the standard FHN equation in order to control the pacemaker period independent of the action potential duration. The multiplication factor in the recovery variable, w, allows us to control the duration of the plateau and the resting state independently by parameters w_H and w_L respectively.

We use D = 0.2 cm² s⁻¹,  = 0.6 s, β = 0.7 and γ = 0.5. Neumann (no-flux) boundary condition is used along the edges of the simulation domain. In the homogeneous excitable medium, the pacemaker current I_P = 0 and w_H = w_L = 0.06 s⁻¹. See figure 1 for the single cell action potential and the phase plane of this model. In a two dimensional homogeneous medium with this set of parameter values, the medium can support stably rotating spiral waves with period ≈0.73 s. The pacemaker is defined as a central disc with radius r_CP in which I_P = 1. We use w_H = 0.06 s⁻¹ and w_L = w_Lcp inside the pacemaker. With all other parameter values fixed, the pacemaker period, T_CP, can be controlled by varying w_Lcp according to the relationship, $T_{\mathrm {CP}}=T_0+\frac {\theta }{w_{\mathrm {Lcp}}}$ where T₀ ≈ 0.3908 s and θ ≈ 0.0460 [28]. The hyperbolic dependence on w_L occurs since w_L controls the rate of depolarization of the action potential. The spiral tip was defined as the site (x*,y*) where v(x*,y*),w(x*,y*) have values close to center of the phase plane (v* = 0,w* = 0) (see figure 1), as demonstrated in Winfree [4]. A two dimensional spline interpolation scheme was used to find the site (x,y) which has the minimum distance to (x*,y*). The interbeat interval (IBI) or the spiral period is the time between successive crossing of v = 0.4 with a positive slope at a given spatial location. In this paper we report IBIs measured from the spiral core (denoted as site i, x = 3.5 mm, y = 5.5 mm) and at a point away from the core (denoted as site ii, x = 5.5 mm, y = 5.5 mm).

Zoom In Zoom Out Reset image size

The equations are integrated using forward Euler discretization with Δt = 0.5 μs on a 200 × 200 grid with Δx = 50 μm. The size of the simulation domain in dimensioned units is 10 × 10 mm². We have crosschecked representative numerical results with simulation using fourth order Runge–Kutta numerical scheme, and further with simulations with Δt = 0.125 μs on a 400 × 400 grid with Δx = 25 μm.

The spiral wave is generated using a broken wave front that was allowed to evolve for 100 s in a medium with w_L = 0.06 s⁻¹. This state is then used as an initial condition for later simulations.

The rotation period of the spiral and its stability are controlled by the parameter w_L (figure 2(a)). Spiral tip trajectories for w_L = 0.036, 0.048 and 0.06 s⁻¹, are illustrated in the insets of figure 2(a). The wave solution is unstable when w_L ⩽ 0.003 s⁻¹. For w_L ⩾ 0.003 s⁻¹ the tip meanders in a flower-like trajectory. As we increase w_L the medium becomes more excitable, and the spiral wave rotates faster. The amplitude of tip movement decreases with increasing w_L, when the system excitability increases [4, 11]. We use w_L = 0.06 in the excitable medium, and thus in the absence of the pacemakers, the tip meanders as shown in the figure 2(a) (see also the supplementary movie 1, available from stacks.iop.org/NJP/15/023028/mmedia). Figure 2(b) shows the IBIs from a point near the spiral tip for w_L = 0.06 s⁻¹. The time series is quasiperiodic locally near the tip but it is periodic with a period 0.73 s away from it [9].

Zoom In Zoom Out Reset image size

In subsequent simulations we introduce a pacemaker in the center with w_L = w_Lcp. When w_Lcp > 0.06 s⁻¹, the pacemaker is faster than the reentrant wave, and it can drive the reentrant spiral out of the medium [30–32]. But when the period of the pacemaker is longer than the period of the reentrant wave, it may still modify the dynamics of the spiral. Let τ denotes the ratio of the spiral period to the period of the pacemaker. We consider 0.012 ⩽ w_L ⩽ 0.048 s⁻¹, corresponding to 0.26 ⩽ τ ⩽ 0.75, and 0.2 < r_CP < 2 mm.

Initially the spiral tip is located at position (x,y) = (3.5,5.5) mm. When the spiral tip is distant from the pacemaker two possible types of behavior can occur. (i) The spiral waves act to reset the pacemaker leading to the entrainment of the pacemaker, but the spiral tip remains distant from the pacemaker. For example, figure 3(a) shows an example of 1:1 entrainment of the pacemaker by the spiral when τ = 0.26 and r_CP = 0.20 mm. (ii) The spiral waves migrate towards the pacemaker and eventually leads to dynamics in which asymptotically the spiral tip is inside or adjacent to the pacemaker boundary for some or all of the time (figure 3(b)). Although these situations are worthy of further investigation, in the current work we focus attention on situations in which the spiral tip has direct interaction with the pacemaker.

Zoom In Zoom Out Reset image size

We consider three types of direct spiral pacemaker interaction:

Type I. The spiral tip is inside the pacemaker and the dynamics depend on the period of the pacemaker. The spiral either rotates rigidly or meanders so that the trajectory of the spiral core is either circular or has a petal geometry.

Type II. The spiral tip stays on the periphery of the pacemaker or is exterior to the pacemaker.

Type III. The spiral tip executes complex trajectories in which it is partially in the interior of the pacemaker and partially on either the circumference or the exterior of the pacemaker.

Regions of different types of dynamics as a function of τ and r_CP are summarized in figure 4(a). This is a schematic representation for which we varied r_CP from 0.2 to 2.0 mm and w_Lcp from 0.012 to 0.048 s⁻¹ in steps of 0.006 s⁻¹ (corresponding to τ values from 0.26, 0.37, 0.47, 0.55, 0.62, 0.69, and 0.75). When r_CP = 0.4 mm, for τ < 0.6, we observe type III behavior, i.e., the tip trajectories move inside and outside the pacemaker. This can lead to quasiperiodic or irregular dynamics in the medium. Type II behavior is observed when τ = 0.26, and r_CP > 0.4 mm (inwardly rotating petals) or when τ = 0.75, with outwardly rotating petals. For intermediate values of τ and r_CP one observes type I dynamics.

Zoom In Zoom Out Reset image size

When the spiral tip meanders, the period of the spiral wave measured in the medium depends on the distance of the measuring site from the spiral tip. For example, in the homogeneous medium the time series measured from a site near a meandering tip is quasiperiodic, but it is nearly periodic if measured from a site away from the tip [9]. We take time series from x = 7.5 mm, y = 7.5 mm (site ii), away from the spiral tip and the pacemaker, to represent the global dynamics in the medium. In the homogeneous medium variation in the spiral period at this location is less than 2% of the mean period. The mean spiral period in the medium is also a function of the pacemaker period (τ) and its radius (r_CP). Figure 4(b) shows that, for τ = 0.26 and 0.37, the mean spiral period from site ii increases linearly with r_CP, similar to what happens during spiral attachment with non-excitable boundary [16]. But for larger τ the spiral period in the medium depends more on τ and less on r_CP. For example, in figure 4(b), when τ is 0.47 and 0.55, the mean spiral period increases with r_CP but then becomes nearly constant for r_CP > 1.4 mm. Then for τ = 0.69 the mean spiral period is independent of r_CP, and for τ = 0.75 the spiral is faster than in the homogeneous medium. The best straight line passing through the linear regime is also shown in figure 4(b). The goodness of fit values (R²) for the linear regime is as follows: 0.998 (τ = 0.26), 0.996 (τ = 0.37), 0.991 (τ = 0.47), 0.984 (τ = 0.55), 0.969 (τ = 0.62) and 0 for (τ = 0.69 and 0.75).

However the fluctuation from the mean period in presence of the pacemaker varies as much as 50%, depending on the type of dynamics (see the supplementary figure S1, available from stacks.iop.org/NJP/15/023028/mmedia for standard deviations from the mean in the distribution of IBIs). Variation in the spiral period is more pronounced for type II and type III dynamics, but it is less than 5% for type I dynamics. Examples of spiral tip trajectories and associated IBIs and IBI return plots are given in figures 4(c)–(h). The IBI and its return plots are calculated from a time series measured at a point away from spiral tip and the pacemaker (at site ii). Details of the dynamics in each of these regions are described below.

3.1. Type I: spiral inside the pacemaker

3.2. Type II: spiral attachment at high τ, Doppler effect

3.3. Type II: spiral attachment at low τ

3.4. Type III: torus doubling and irregular dynamics

3.5. Multistability

The effect of spiral wave interactions with unexcitable boundaries decay exponentially as a function of the distance of the spiral tip from the boundary [34]. Similarly the effect of the pacemaker on the wave dynamics depends on the relative position of the spiral tip with respect to the pacemaker. For example, consider the initial condition employed in this study, figure 5(a). With τ = 0.75 and r_CP = 1.8 mm, the final dynamics for this initial condition is type II. However, if the spiral tip was closer to the pacemaker the tip can enter inside the pacemaker, resulting in type I dynamics. On the other hand, if the initial location of the spiral tip was farther, the spiral wave is unaffected by the pacemaker.

Zoom In Zoom Out Reset image size

To illustrate the initial condition dependence on the spiral pacemaker interaction, we show dynamics starting from two different initial conditions. Consider a pacemaker with τ = 0.27 and r_CP = 1.8 mm located at (1,1) mm. The initial position of the spiral as in figure 5(a) is the initial condition used throughout this study. Let this initial condition be IC₀. Starting from IC₀ with the specified pacemaker will lead to type II dynamics, and the final tip trajectory is as shown in the first panel of figure 5(e). Call this state IC₁. Using the spiral position as in IC₁ as the initial condition, we vary τ and let the dynamics continue for 10 s. The resulting tip trajectories for τ =  0.37, 0.46, 0.55, 0.62, 0.69 and 0.75 are shown in figure 5(e). Similarly, starting from IC₀ with the pacemaker of same size (r_CP = 1.8 mm) but with τ = 0.75, the final dynamics is type I, with a tip trajectory as shown in the last panel of figure 5(f). Let this be IC₂. The tip trajectories obtained by using IC₂ as the initial condition for pacemakers with τ varying from, 0.27 to 0.69 are shown in the first six panels of figure 5(f). In this example the final tip trajectories are independent of the given initial conditions (i.e. IC₀, IC₁ and IC₂) for intermediate values of τ (0.26 < τ < 0.75). But for τ = 0.27, the final dynamics is type II if started from IC₀, and type I if started from IC₂. Similarly, for τ = 0.75, the final dynamics is type I if started from IC₀, and type II if started from IC₁.

Our study shows that local heterogeneities like pacemakers can change spiral wave dynamics in an excitable medium and can lead to quasiperiodic or irregular dynamics. In particular, when the pacemaker is slower than the spiral wave three different types of dynamics can occur, depending on the radius and period of the pacemaker. The spiral wave can move inside the pacemaker (type I), can drift and meander along the boundary of the pacemaker (type II), or it can move partly inside and partly outside the pacemaker (type III). Type I dynamics can lead to periodic or quasiperiodic activation, depending on the period of the pacemaker. In type II dynamics, activation in the medium is quasiperiodic while in type III, the dynamics can also be irregular. The types of dynamics as a function of period and the radius of the pacemaker are summarized in the schematic diagram (figure 4(a)).

There have been a large number of earlier studies concerned with the interaction of spiral waves with heterogeneities. It is well known that spiral waves can become anchored to localized heterogeneities [15, 16, 18, 19], and that global gradients of excitability in the medium can lead to the drifting of spiral waves [35]. The work most closely related to ours involves the study of the effect of localized partially excitable heterogeneities on spiral wave dynamics. As in the current work, localized heterogeneities can lead to the drift of spiral tips. Asymptotically, the spiral tip can be entirely internal to the heterogeneity or remain on the periphery [20, 22–25], similar to what we have observed in figure 4.

In addition to these behaviors, in the current work we have found that the spiral tip can display complex dynamics that can be partially characterized by considering the return map showing the dependence of the IBI in the excitable medium on the preceding IBI (figure 4). Unlike the situation during quasiperiodic dynamics when the return map is a closed curve, we have found two variants. For R = 0.47 r_CP = 0.4 mm, figure 4(g), there appears to be a structure similar to those observed in period doubling bifurcations on a torus [33] and also during spiral wave activity in the Luo–Rudy model of cardiac wave propagation [9]. When R = 0.26 r_CP = 0.4 mm, the dynamics are more irregular and the IBI return map appears random (figure 4(h)). Although this is qualitatively similar to what was observed during spiral breakup [9], in the current situation the spiral wave does not break up.

These observations may be relevant to clinically observed cardiac arrhythmias. The cardiac rhythm is normally set by a localized pacemaker called the sinus node located in the right atrium. In sinus node reentrant tachycardia there is a reentrant wave centered in the vicinity of the sinus node. It is a relatively rare arrhythmia and is difficult to characterize [36, 37]. Further work is needed to clarify if rhythms similar to those described here are observed in sinus node reentrant tachycardia or other reentrant tachycardias that appear to circulate on a micro rather than a macro reentrant pathway.

Independent of the clinical utility, the current paper highlights a range of complex dynamics that are expected in excitable media with embedded localized pacemakers. The mathematical model used simplifies underlying structural and biophysical details of cardiac pacemakers. It will be interesting to compare our results with those from models which incorporate anatomical details and ionic mechanisms. This would help identify both basic and detailed ingredients important for spiral wave generation, maintenance and termination. It will also be of interest to better characterize the boundaries for the different dynamical regimes and to compute the bifurcations that arise as the size and period of the pacemaker are varied.

We thank NSERC and the Centre for Applied Mathematics in Bioscience and Medicine (McGill University) for financial support.