## Abstract

Heterogeneity of cardiac tissue is an important factor determining the initiation and dynamics of cardiac arrhythmias. In this paper, we studied the effects of gradients of electrophysiological heterogeneity on reentrant excitation patterns using computer simulations. We investigated the dynamics of spiral waves in a two-dimensional sheet of cardiac tissue described by the Luo-Rudy phase 1 (LR1) ventricular action potential model. A gradient of action potential duration (APD) was imposed by gradually varying the local current density of K^{+} current or inward rectifying K^{+} current along one axis of the tissue sheet. We show that a gradient of APD resulted in spiral wave drift. This drift consisted of two components. The longitudinal (along the gradient) component was always directed toward regions of longer spiral wave period. The transverse (perpendicular to the gradient) component had a direction dependent on the direction of rotation of the spiral wave. We estimated the velocity of the drift as a function of the magnitude of the gradient and discuss its implications.

- ionic model
- spiral wave
- reentrant arrhythmias
- tissue heterogeneity

various experimental and theoretical investigations show that many dangerous cardiac arrhythmias are driven by reentrant sources of excitation (1, 10, 29, 30). In many of these cases, such reentrant sources are shown to be rotating spiral waves of excitation (5, 6, 11, 16, 21, 26). The dynamics of spiral waves are considered to be an important factors in determining the type of cardiac arrhythmia that will occur: stationary rotation of spiral waves is associated with monomorphic tachycardias, whereas nonstationary rotation or drift can cause polymorphic tachycardias and torsade de pointes (9).

Spiral wave dynamics are in many ways affected by the ionic heterogeneity of cardiac tissue. Ionic heterogeneity is considered to be one of the main factors underlying the initiation of spiral waves in the heart (15, 19). In addition, ionic heterogeneity is thought to play an important role in the transition from a spiral wave excitation pattern to the spatiotemporally irregular pattern characteristic of fibrillation (13, 19).

There are many types of ionic heterogeneity in the heart. Well-known examples are the base-apex (3) and endocardial-epicardial (3, 17, 32) action potential duration (APD) gradients in the ventricles and differences in APD between left and right ventricles (24). Other examples are the APD gradient running from the crista terminalis to pectinate muscles in the right atrium (25) and the difference in APD and effective refractory period between left and right atria (25).

Ionic heterogeneity in the heart increases under particular conditions, such as ischemia or chronic heart failure, or in genetic disorders, such as those associated with long QT and Brugada syndromes (2, 27). In most of these cases, the increased ionic heterogeneity is associated with an increased occurrence of cardiac arrhythmias.

There have been several studies on the effects of stepwise ionic heterogeneities of cardiac tissue on spiral wave dynamics. Xie et al. (31) showed in a model study that stepwise heterogeneities could induce spiral breakup even if the action potential restitution slope is shallow. Experimental studies conducted by Fast and Pertsov (7) showed that spiral waves drift along the border of a stepwise ionic heterogeneity induced by the local application of quinidine.

Another important type of organization of ionic heterogeneity is in the form of a smooth gradient. Such heterogeneities can occur either as a result of smooth variation of ionic properties of cardiac tissue or as a result of the smoothing effect of electrotonic interactions between cardiac cells, even if properties change abruptly from one cell to the next (14, 28). A well-known example of the latter type of heterogeneity is the endocardial-epicardial gradient across the ventricular wall, which is largely caused by the presence of the electrophysiologically distinct M cells (17).

So far, only research using simplified FitzHugh-Nagumo models of excitable tissue has been done on the effects of smooth ionic gradients on the dynamics of spiral waves (22). In these simplified models, spiral waves had circular cores and action potential shapes very different from those occurring in cardiac tissue. Therefore, the purpose of the present paper was to study the effect of gradients of electrophysiological properties on the dynamic behavior of spiral waves in computer simulations of a realistic ionic model of cardiac tissue. To model cardiac tissue, we chose the Luo-Rudy phase 1 (LR1) model for ventricular cardiomyocytes (18), which was used previously by Xie et al. (31) to study the effects of stepwise ionic heterogeneities.

We show that, as in the case of FitzHugh-Nagumo models, drift can be decomposed into a longitudinal component, parallel to the gradient, and a transverse component, perpendicular to the gradient. Our main finding is that independent of the type of electrophysiological gradient, the longitudinal component of drift is always directed toward regions of longer spiral wave period. We found drift velocities on the order of 0.002 mm/ms for period gradients of 0.2 ms/mm, amounting to a displacement of 2 mm during 1 s of spiral rotation.

These results may help in both the predicting and understanding of the behavior of reentrant wave patterns in the proximity of ischemic zones and other regions with a well-known effect on spiral period.

## MATERIALS AND METHODS

#### Mathematical modeling.

We used an ionic model to study the propagation of waves of excitation in cardiac tissue. Ignoring the discrete character of microscopic cardiac cell structure, cardiac tissue can be modeled as a continuous system using the following partial differential equation (PDE)
Equation 1where *V* is membrane potential (in mV), *t*is time (in ms), *C*
_{m} = 1 μF/cm^{2} and is membrane capacitance, *D* = 0.001 cm^{2}/ms and is the diffusion coefficient, and*I*
_{ion} is the sum of all transmembrane ionic currents (in mA/cm^{2}). For *I*
_{ion}, we used the following equation: *I*
_{ion} =*I*
_{Na} + *I*
_{si} +*I*
_{K} + *I*
_{K1} +*I*
_{Kp} + *I*
_{b} as described in the LR1 model (18), where*I*
_{Na} is Na^{+} current,*I*
_{si} is slow inward current,*I*
_{K} is K^{+} current,*I*
_{K1} is inward rectifying K^{+} current,*I*
_{Kp} is K^{+} pump current, and*I*
_{b} is background current. Here,*I*
_{Na} =*G*
_{Na}
*m*
^{3}
*hj*(*V*− *E*
_{Na}), *I*
_{si} =*G*
_{si}
*df*(*V* −*E*
_{si}), *I*
_{K} =*G*
_{K}
*xx*
_{1}(*V*− *E*
_{K}), *I*
_{K1} =*G*
_{K1}
*K*
_{1∞}(*V*− *E*
_{K}), *I*
_{Kp} =*G*
_{Kp}(*V* −*E*
_{K}), and *I*
_{b} =*G*
_{b}(*V* −*E*
_{b}), where *G* is conductance,*E* is Nernst potential, and *K*
_{100}models the inward rectification of *I*
_{K1}. The gating variables *m*, *h*, *j*,*d*, *f*, and *x* are governed by a Hodgkin-Huxley-type differential equation.

The parameter settings were as in the original LR1 model except for the*G*
_{si}, *G*
_{K}, and*G*
_{K1} conductances. For*G*
_{si}, values of 0, 0.030, 0.035, 0.040, and 0.045 were used to vary the meander pattern (31). The upper value of 0.045 was chosen to avoid spiral breakup occurring. Under homogenous tissue conditions, *G*
_{K} was set to 0.705 to shorten APD (31) and *G*
_{K1}was set to 0.6047 as in the original LR1 model.

We studied the effects of gradients of electrophysiological properties caused by local differences in *I*
_{K} and*I*
_{K1} densities. To simulate heterogeneity of*I*
_{K1} density, a *G*
_{K1}gradient was applied by varying *G*
_{K1} linearly from a value of 0.423 to a value of 0.605. Similarly, to mimic heterogeneity of *I*
_{K} density,*G*
_{K} was varied linearly from a value of 0.600 to a value of 0.705. Gradients were applied such that the parameter value increased in a positive *y*-direction. To induce gradients with different slopes, tissue sizes were varied from 250 × 250 (5 × 5 cm) to 700 × 700 (14 × 14 cm) nodes.

#### Computer simulations.

Two-dimensional tissue was simulated by integrating the PDE described in *Eq. 1
*.

To speed up computations, reaction and diffusion were separated using operator splitting. The diffusion PDE was solved using a time step of Δ*t* = 0.1 ms. The reaction ordinary differential equation was solved using a time-adaptive forward Euler scheme with two different time steps: Δ*t*
_{large} = 0.1 ms and Δ*t*
_{small} = 0.02 ms. By default, the ordinary differential equation was solved using Δ*t*
_{large}. However, if ∂*V*/∂*t* > 1, the results were discarded and computations were repeated iterating five times over Δ*t*
_{small}. The equations for the gating variables were integrated using the Rush and Larsen scheme (23). In all simulations, the space step was set to Δ*x* = 0.02 cm. We checked the accuracy of our variable time step integration method by comparing it with the conventional Euler integration scheme for a selected subset of the simulations and found similar results (data not shown).

Spiral wave reentry was initiated by applying a S1-S2 protocol with parallel electrode positioning. To establish meander patterns and drift, spiral tip trajectories were traced using the algorithm presented by Fenton and Karma (8) along an isopotential line of *V* = −35 mV.

All simulations were written in C++ and run on an Intel Pentium III personal computer with an 800-MHz central processing unit.

## RESULTS

#### Spiral dynamics in homogeneous tissue.

Figure 1
*A* shows the typical dynamics of a spiral wave in the LR1 model. The tip of the spiral wave followed a meandering trajectory with a distinctive hypocycloidal pattern made up of five outward petals.

Changing the excitability of cardiac tissue by changing*G*
_{si} resulted in different spiral dynamics, as shown in Fig. 2. For*G*
_{si} = 0 (low excitability), the spiral core was circular; for *G*
_{si} = 0.030, the core had a hypocycloidal pattern similar to the one shown in Fig. 1. For ever higher excitability (from *G*
_{si} = 0.035 to*G*
_{si} = 0.045), this hypocycloidal pattern persisted; however, its size increased with increasing*G*
_{si} values.

#### Spiral dynamics in tissue with a gradient of heterogeneity.

We studied the dynamics of spiral waves in tissue with a gradient of heterogeneity. A gradient was created by gradually varying the local current density of either *I*
_{K} or*I*
_{K1} by varying *G*
_{K} and*G*
_{K1} values, respectively, similar to the approach taken by Xie et al. (31).

Figure 3 shows a series of snapshots of spiral wave dynamics in tissue with a *G*
_{K}gradient ranging from *G*
_{K} = 0.600 at the bottom of the medium to *G*
_{K} = 0.705 at the top of the medium. One can see that, over the course of time, the spiral wave gradually shifted from its initial position (indicated by the white quadrant) to a new position located down and to the right of that initial position.

Figure 4
*A* shows the tip trajectory of the spiral from Fig. 3. Because of substantial meandering, the tip trajectory looks rather complicated; however, the net downward and to the right drift can be clearly seen. To separate the drift from the meandering, we averaged the trajectory from Fig.4
*A* over the characteristic time of the meandering pattern (300 ms). We saw (Fig. 4
*B*) that the spiral drift occurred at an approximately constant angle relative to the *y*-axis. Therefore, spiral drift can be regarded as a vector consisting of two components: one parallel to the gradient of heterogeneity (along the*y*-axis), and the other perpendicular to this gradient. We call these the longitudinal and transverse components of drift, respectively.

In Figs. 3 and 4, the value of *G*
_{K} increases with increasing *y*-coordinates, i.e., APD and hence the spiral wave period decreases along the positive *y*-direction. It thus follows that the longitudinal component of spiral drift was directed toward regions of longer period (Fig.4
*C*).

The direction of the transverse component of drift depends on the direction of rotation of the spiral wave (22). In Fig.4
*D*, we show drift of a spiral with an opposite direction of rotation under the same parameter settings as in Fig. 4
*B*. The longitudinal component of drift did not changed, whereas the transverse component of drift had an opposite direction. From a mathematical point of view, this result is trivial as it follows from the symmetry of the system. The direction of transverse drift can be written as
Equation 2where **
** and **
** are unit vectors along the transverse and longitudinal drift components, respectively, and **
** is the unit vector of the angular velocity of the spiral wave (22). From *Eq. 2
*, it follows that if the direction of rotation reverses (**
**), the direction of transverse drift (**
**) also reverses.

Similar computations were done with a gradient of*I*
_{K1} density. All other parameters were kept the same as in Fig. 4. One can see (Fig. 5) that despite the different gradient, qualitative characteristics of drift remained the same: drift had two components, a longitudinal component, directed toward regions with longer spiral wave period; and a transverse component, with its direction given by *Eq. 2
*.

In Figs. 4 and 5, the spiral wave had a meander pattern corresponding to a *G*
_{si} value of 0.030. To investigate how the meander pattern affects spiral wave drift, we studied drift for all five meander patterns shown in Fig. 2. In addition, we investigated the influence of gradient magnitude on drift of the spiral wave.

Figure 6 shows the speed of the longitudinal drift of spiral waves for different meander patterns (values of *G*
_{si}) as a function of*G*
_{K} or *G*
_{K1} gradient magnitudes. We saw that independent of the type (*G*
_{K} or *G*
_{K1}) and magnitude of the gradient and independent of the meander pattern of the spiral wave, the longitudinal drift was always directed toward regions of longer spiral wave period. It can also be seen that velocity increased approximately linearly with increasing magnitude of the gradient. In addition, we see in Figs. 4 and 5 that for higher values of*G*
_{si}, longitudinal drift speed increases and the graphs have a steeper slope.

Figure 7 shows the dependence of transverse drift speed on the meander pattern and gradient magnitude. One can see that for most values of *G*
_{si}, transverse drift velocity was positive, i.e., it coincided with the direction it had in Fig. 4
*B*. However, for the lowest excitability (*G*
_{si} = 0), transverse speed was virtually absent for the *G*
_{K1} gradient and was very small and even had an opposite direction for the*G*
_{K} gradient.

For intermediate values of *G*
_{si} (from 0.030 to 0.040), transverse drift speed showed the same linear dependence on gradient magnitude as longitudinal drift speed. However, for the lowest (*G*
_{si} = 0) and highest excitability (*G*
_{si} = 0.045), this dependency was substantially nonlinear. Also, the influence of the*G*
_{si} value on transverse drift speed was less strong than for longitudinal drift speed and displayed saturation from*G*
_{si} = 0.040 onward.

From a comparison of Figs. 6 and 7, it follows that transverse drift speed was a factor two to three times smaller than longitudinal drift speed. In addition (looking only at intermediate values of*G*
_{si}), it follows that the slope of the linear dependence of drift speed on magnitude of the gradient was less steep for transverse than for longitudinal drift. This implies that the angle of drift relative to the gradient decreased for increasing gradient magnitude.

In all cases studied above, the spiral waves drifted to regions of longer spiral wave period and did so faster for stronger gradients. This suggests that the gradient in period (induced by a*G*
_{K} or *G*
_{K1} gradient) was the driving force behind the drift. In addition, we saw that the drift velocity also depended on the meandering pattern, which was varied by changing the parameter *G*
_{si}. However,*G*
_{si} not only influences the meander pattern (Fig. 1) but also influences heterogeneity: for different values of*G*
_{si}, the same gradients in*G*
_{K} or *G*
_{K1} produce different gradients in spiral wave period. This raises the question of whether *G*
_{si} affects the spiral drift by changing the gradient or by changing the dynamics of spiral wave meandering.

To answer the above question, we plotted the longitudinal drift speeds from Fig. 6 and the transverse drift speeds from Fig. 7 as a function of the magnitude of the effective period gradient induced by a particular combination of the *G*
_{K} or*G*
_{K1} gradient magnitude and*G*
_{si} value. The resulting graphs are shown in Fig. 8. We can clearly see that the longitudinal drift velocities for the different values of G_{si} and the different gradient types have converged toward one another: the graphs have a similar slope and similar (extrapolated) intersection point with the *y*-axis. This suggests that longitudinal drift speed was linearly correlated to period gradient magnitude and that neither the meander pattern nor type of gradient played a role here.

For transverse drift, the situation is less obvious. For intermediate values of *G*
_{si}, transverse drift speed graphs also converged, suggesting a linear dependence on period gradient magnitude. However, for the lowest (*G*
_{si} = 0) and highest excitability (*G*
_{si} = 0.45), the graphs for both types of gradient lie below that for the intermediate excitability levels (from *G*
_{si}= 0.030 to *G*
_{si} = 0.040). These results might indicate that transverse drift speed was not determined by period gradient magnitude alone. The meander pattern could have also played an important role here.

Note that all results were checked for dependence on initial conditions and timing of S1-S2 stimulation. Found differences in longitudinal and transverse drift speeds stayed within 1%, indicating that drift behavior did not depend on the initial phase of rotation.

#### Figure-eight reentry in a gradient of heterogeneity.

A frequently observed reentrant pattern is so-called figure-eight reentry (4), which is composed of two interconnected, counterrotating spiral waves. Because the direction of transverse drift in a gradient of heterogeneity depends on the direction of spiral wave rotation relative to the gradient direction, two counterrotating spirals should either diverge or converge (12, 20). Figure9 shows drift patterns of figure-eight reentries in a medium with a *G*
_{K1} gradient.

With the use of a S1-S2 protocol with different S1 electrode positions, we generated two differently oriented figure-eight reentries. As predicted by the results above, in Fig. 9, *A–C*, the two spirals diverged, whereas in Fig. 9, *D–F*, they converged. In both cases, the figure-eight reentry moved downward, i.e., toward a longer spiral period. Note that convergence of spirals did not lead to mutual annihilation. Instead, we observed that after spirals had converged to some minimal distance, no further convergence occurred and only the longitudinal component of drift remained.

## DISCUSSION

We studied the effect of gradients of heterogeneity on the dynamics of spiral waves in the LR1 model. It was demonstrated that spiral waves drift in the presence of a gradient of heterogeneity and that this drift consists of two components: a longitudinal component (parallel to the gradient) and a transverse component (perpendicular to the gradient).

Longitudinal drift was always directed toward regions of longer spiral wave period. This result is similar to findings obtained using a two-equation FitzHugh-Nagumo model (22). The fact that this result generalizes across such very different models enables us to assume that it is a general phenomenon that should exist in other models of cardiac tissue as well as in experiments. In addition, we found that the velocity of longitudinal drift was linearly proportionate to the magnitude of the spiral period gradient and independent of the spiral wave meander pattern.

The direction of the transverse component of drift was shown to be given by *Eq. 2
*. This finding agrees with data from Fast and Pertsov (7), who studied spiral wave dynamics under a stepwise heterogeneity in an experimental setup. For intermediate values of excitability, transverse speed was also linearly proportionate to the period gradient magnitude. For the lowest and highest excitability cases, the dependency was rather different, suggesting that the meander pattern might play a role here.

One of the interesting implications of our findings is that we can predict how the dynamics of a spiral wave will be influenced by a particular gradient present in the heart. The ventricular base-apex gradient, for example, should cause a drift of transmurally oriented spirals toward the septum, where APD is longest (25). The transmural endocardial-epicardial gradient should cause a drift of intramurally oriented spirals toward the midmyocardial region, where APD is longest due to the presence of M cells (17). Note, however, that the influence of other factors, such as the three-dimensional nature of reentry in the ventricles and the presence of rotational anisotropy, should also be taken into account. Similar predictions can be made for spiral behavior in the proximity of an ischemic border zone, where a spiral wave should move away from the ischemic region, where APD is shortest. Note, however, that because multiple electrophysiological factors change during ischemia, spiral behavior under these conditions requires further study.

The drift velocity found in our computations is in the order of 0.002 mm/ms in a period gradient of 0.2 ms/mm. During a single second of drift, a spiral wave thus can travel a distance of ∼2 mm, which can have a significant effect on cardiac arrhythmias and the appearance they make on an ECG.

In conclusion, the aim of the present study was to investigate the basic effects of gradients of ionic heterogeneity in cardiac tissue on spiral wave dynamics. We showed that these gradients lead to drift of spiral waves toward regions of longer period independent of the type of ionic heterogeneity. Our main conclusion is that differences in spiral wave period are the driving force behind the drift of spiral waves.

## Acknowledgments

This research was supported by Netherlands Organization for Scientific Research Grant 620061351 of the Research Council for Physical Sciences.

## Footnotes

Address for reprint requests and other correspondence: K. H. W. J. ten Tusscher, Dept. of Theoretical Biology, Utrecht Univ., Padualaan 8, 3584 CH Utrecht, The Netherlands (E-mail: khwjtuss{at}hotmail.com).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.First published October 10, 2002;10.1152/ajpheart.00608.2002

- Copyright © 2003 the American Physiological Society