Anisotropic shortening in the wavelength of electrical waves promotes onset of electrical turbulence in cardiac tissue: An in silico study.

Several pathological conditions introduce spatial variations in the electrical properties of cardiac tissue. These variations occur as localized or distributed gradients in ion-channel functionality over extended tissue media. Electrical waves, propagating through such affected tissue, demonstrate distortions, depending on the nature of the ionic gradient in the diseased substrate. If the degree of distortion is large, reentrant activity may develop, in the form of rotating spiral (2d) and scroll (3d) waves of electrical activity. These reentrant waves are associated with the occurrence of lethal cardiac rhythm disorders, known as arrhythmias, such as ventricular tachycardia (VT) and ventricular fibrillation (VF), which are believed to be common precursors of sudden cardiac arrest. By using state-of-the-art mathematical models for generic, and ionically-realistic (human) cardiac tissue, we study the detrimental effects of these ionic gradients on electrical wave propagation. We propose a possible mechanism for the development of instabilities in reentrant wave patterns, in the presence of ionic gradients in cardiac tissue, which may explain how one type of arrhythmia (VT) can degenerate into another (VF). Our proposed mechanism entails anisotropic reduction in the wavelength of the excitation waves because of anisotropic variation in its electrical properties, in particular the action potential duration (APD). We find that the variation in the APD, which we induce by varying ion-channel conductances, imposes a spatial variation in the spiral- or scroll-wave frequency ω. Such gradients in ω induce anisotropic shortening of wavelength of the spiral or scroll arms and eventually leads to instabilitites.

Introduction

Nonlinear waves in the form of spirals occur in many excitable media, examples of which include Belousov-Zhabotinsky-type systems [1], calcium-ion waves in Xenopus oocytes [2], the aggregation of Dictyostelium discoideum by cyclic-AMP signaling [3], the oxidation of carbon monoxide on a platinum surface [4], and, most important of all, cardiac tissue [5]. Understanding the development of such spiral waves and their spatiotemporal evolution is an important challenge in the study of extended dynamical systems, in general, and especially in cardiac tissue, where these waves are associated with abnormal rhythm disorders, which are also called arrhythmias. Cardiac tissue can support many patterns of nonlinear waves of electrical activation, like traveling waves, target waves, and spiral and scroll waves [6]. The occurrence of spiral- and scroll-wave turbulence of electrical activation in cardiac tissue has been implicated in the precipitation of life-threatening cardiac arrhythmias like ventricular tachycardia (VT) and ventricular fibrillation (VF), which destroy the regular rhythm of a mammalian heart and render it incapable of pumping blood. These arrhythmias are the leading cause of death in the industrialized world [7-11].

Biologically, VF can arise because of many complex mechanisms. Some of these are associated with the development of instability-induced spiral- or scroll-wave turbulence [12]. One such instability-inducing factor is ionic heterogeneity [13, 14], which arises from variations in the electrophysiological properties of cardiac cells (myocytes), like the morphology and duration of their action-potentials (APs) [15-18]. Such variations may appear in cardiac tissue because of electrical remodeling [19-21], induced by alterations in ion-channel expression and activity, which arise, in turn, from diseases [22] like ischemia [23, 24], some forms of cardiomyopathy [25], and the long-QT syndrome [26]. To a certain extent, some heterogeneity is normal in healthy hearts; and it has an underlying physiological purpose [16, 27–31]; but, if the degree of heterogeneity is more than is physiologically normal, it can be arrhythmogenic [20, 24, 32]. It is important, therefore, to explore ionic-heterogeneity-induced spiral- or scroll-wave turbulence in mathematical models of cardiac tissue, which allow us to control this heterogeneity precisely, in order to be able to identify the nonlinear-wave instability that leads to such turbulence. We initiate such a study by examining the effects of this type of heterogeneity in three cardiac-tissue models, which are, in order of increasing complexity and biological realism, (a) the two-variable Aliev-Panfilov model [33], (b) the ionically realistic O’Hara-Rudy (ORd) model [34] in two dimensions (2D), and (c) the ORd model in an anatomically realistic simulation domain. In each one of these models, we control parameters (see below) in such a way that the ion-channel properties change anisotropically in our simulation domains, thereby inducing an anisotropic spatial variation in the local action potential duration APD. We show that this variation in the APD leads, in all these models, to an anisotropic reduction of the wavelength of the spiral or scroll waves; and this anisotropic reduction of the wavelength paves the way for an instability that precipitates turbulence, the mathematical analog of VF, in these models.

Materials and methods

The Aliev-Panfilov model provides a simplified description of an excitable cardiac cell [33]. It comprises a set of coupled ordinary differential equations (ODEs), for the normalized representations of the transmembrane potential V and the generalized conductance r of the slow, repolarizing current: fast processes are governed by the first term in Eq (1), whereas, the slow, recovery phase of the AP is determined by the function in Eq (2). The parameter a represents the threshold of activation and k controls the magnitude of the transmembrane current. We use the standard values for all parameters [33], except for the parameter k. We write k = g × k, where g is a multiplication factor and k is the control value of k. In 2D simulations we introduce a spatial gradient (a linear variation) in the value of k along the vertical direction of the domain. To mimic the electrophysiology of a human ventricular cell, we perform similar studies using a slightly modified version of the ionically-realistic O’Hara-Rudy model (ORd) [34, 35]. Here, the transmembrane potential V is governed by the ODE where I, the membrane ionic current, for a generic ion channel x, of a cardiac cell, is where C = 1 μF is the membrane capacitance, f1(p) and f2(p) are, respectively, functions of probabilities of activation (p) and inactivation (p) of the ion channel x, and E is its Nernst potential. We give a list of all the ionic currents in the ORd model in Table 1. We write G = g × G, where G is the original value of the maximal conductance of the ion channel x in the ORd model, and g is a multiplication factor. We model gradients in G as follows: here, L is the length of the side of the square simulation domain, and g and g are, respectively, the maximal and minimal values of g; we can impose gradients in k in the Aliev-Panfilov model in the same manner. For simplicity, we induce the gradient along one spatial direction only: the vertical axis in 2D; and the apico-basal (apex-to-base) direction in 3D. The spatiotemporal evolution of V in both models is governed by the following reaction-diffusion equation: where is the diffusion tensor, and and kV(V − a)(V − 1) + Vr for ORd and Aliev-Panfilov models, respectively. For the numerical implementation of the diffusion term in Eq (6), we follow Refs. [35, 36]. We construct our anatomically realistic simulation domain with processed human-ventricular data, obtained by using Diffusion Tensor Magnetic Resonance Imaging (DTMRI) [37]. For our 2D isotropic-domain with ORd model, we set D = 0.0012cm2/ms. The temporal and spatial resolutions are set to be δx = 0.02 cm and δt = 0.02 ms, respectively, and all the simulations are performed in a domain with 960 × 960 grid points. For the anatomically-realistic domain, we use a phase-field method for the boundary conditions [38]. The value of diffusion constant along the fiber (D∥) is set equal to the value of D in the 2D isotropic case (i.e., 0.0012 cm2/ms) and its value perpendicular (D⊥) to the fiber is 1/4 times D∥. The simulation is performed in a cubical domain with 5123 grid points with the same spatial and temporal resolutions that we use in our 2D simulations. We do not incorporate the intrinsic ionic heterogeneities that are present in real mammalian hearts [16, 27–31]. In our single-cell simulations, the APD is calculated by measuring the duration over which the cell depolarizes and repolarizes to 90% (APD90) of its peak transmembrane voltage in the action potential.

Table 1

The various ionic currents incorporated in the ORd model are tabulated above.

The symbols used for the currents follow Ref. [34].

INa	fast inward Na+ current
Ito	transient outward K+ current
ICaL	L-type Ca2+ current
IKr	rapid delayed rectifier K+ current
IKs	slow delayed rectifier K+ current
IK1	inward rectifier K+ current
INaCa	Na+/Ca2+ exchange current
INaK	Na+/K+ ATPase current
INab	Na+ background current
ICab	Ca2+ background current
IpCa	sarcolemmal Ca2+ pump current
IKb	K+ background current
ICaNa	Na+ current through the L-type Ca2+ channel
ICaK	K+ current through the L-type Ca2+ channel

Results

Spiral-wave instability

In Fig 1(a) we show the variation, with the parameter g, of , where APD is the control APD value for g = 1. We find that decreases with increasing g. Changes in the APD at the single-cell level influence electrical-wave dynamics at the tissue level. In particular, such changes affect the rotation frequency ω of reentrant activity (spiral waves). If θ and λ denote, respectively, the conduction velocity and wavelength of a plane electrical wave in tissue, then , λ ≃ θ × APD. Therefore, if we neglect the effects of curvature [39] and excitable gap, the spiral-wave frequency We find, in agreement with this simple, analytical estimate, that ω decreases as the APD increases. We show this in Fig 1(b) by plotting versus g; here, ω0 is the frequency for g = 1. For the parameter a this simple relation between ω and APD is not observed, because change in a affects not only the APD but also other quantities like θ, which has effects on the value of ω. The spiral-wave frequency ω is obtained by simulating a spiral-wave in a homogeneous domain for every value of g.

Fig 1

Occurrence of spiral-wave instability in a generic model of cardiac tissue with ionic heterogeneity.

Variation of (a) and (b) (see text) with k. ; here, APD is the control value of APD at g = 1 (so at g = 1); we also use other combinations of (, g) in our numerical simulations. We find that decreases with increasing k; however, increases with increasing k. (c) Pseudocolor plots of V, at three representative times (time increases from left to right), illustrating the development of the spiral-wave instability in the Aliev-Panfilov model with a linear gradient in k; S1 Video in SI shows the complete spatiotemporal evolution of this instability.

Similarly, in the ionically realistic ORd model, changes in the ion-channel conductances G alter the APD of the cell and, therefore, the spiral-wave frequency ω. In Fig 2(a1) and 2(a2) we present a family of plots to illustrate the variation in with changes in G. We find that decreases with an increase in g for most currents (I, I, I, I and I); but it increases for some other currents (I, I and I). The rate of change of is most significant when we change G; by contrast, it is most insensitive to changes in G and G. In Fig 2(b1) and 2(b2) we show the variation of with g for different ion channels x. We find that changes in G, which increase APD, decrease ω and vice versa; this follows from Eq (2). The sensitivity of ω, with respect to changes in G, is most for G = G and least for G = G: increases by , as g goes from 0.2 to 5; for G, the same variation in g decreases the value of by . We have done many simulations for each G with different values of to check if a critical value exists such that above (below) we see wave breaks (no wave breaks) for all Gs. We find, however, that no such exists, that is common for all Gs, which is because the stability of the spiral waves depends on the local values of the gradients in APD.

Fig 2

Variation of and with various ionic cunductances.

Plots of and versus g; here, and , where APD = 250 ms, and ω0 = 4.38 Hz are, respectively, the control values of APD and ω; (a1) and (a2) show, respectively, that decreases with the conductances G, for the currents I, I, I, I and I; however, it increases with increasing G, for the currents I, I and I; (b1) and (b2) show that the variation of , with the various channel conductances, is consistent with Eq (2).

We now investigate the effects, on spiral-wave dynamics, of spatial gradients in k, in the 2D Aliev-Panfilov model, and in G, in the 2D ORd model. A linear gradient in k, in the Aliev-Panfilov model, induces a gradient in (see Fig 1(b)); and such a spatial gradient in induces a spiral-wave instability in the low- region. In Fig 1(c) we demonstrate how a gradient in k (g = 1.5 and g = 0.5) leads to the precipitation of this instability (also see S1 Video).

Similarly, for each current listed in Table 1 for the ORd model, we find wave breaks in a medium with a gradient in G. We illustrate, in Fig 3, such wave breaks in our 2D simulation domain, with a gradient (∇G) in any G, for 3 representative currents; we select I, because it has the maximal impact on the single-cell APD, and also on ω in tissue simulations; and we choose I and I, because they have moderate and contrary effects on APD and ω (Fig 2). Our results indicate that gradient-induced wave breaks are generic, insofar as they occur in both the simple two-variable (Aliev-Panfilov) and the ionically realistic (ORd) models of cardiac tissue. In Fig 3(d)–(3f), we present power spectra of the time series of V, recorded from a representative point of the simulation domain; these spectra show broad-band backgrounds, which are signatures of chaos, for the gradients ∇G and ∇G; however, the gradient ∇G induces wave breaks while preserving the periodicity of the resultant, reentrant electrical activity, at least at the points from which we have recorded V.

Fig 3

Spiral-wave instability.

Pseudocolor plots of the transmembrane potential V illustrating spiral-wave instabilities from our numerical simulations of the 2D ORd model for human ventricular tissue, with spatial gradients in (a) G, (b) G, and (c) G (because G decreases with g (Fig 2), the gradient in G must be chosen to be the negative of that in Eq (5)); in (a)-(c) the local value of ω decreases from the top of the simulation domain to its bottom. Power spectra of the time series of V, from representative points in our simulation domain, are shown for gradients in (d) G, (e) G, and (f) G; the spectra in (d) and (e) are consistent with the onset of spiral-wave turbulence; the power spectrum in (f) shows the continuation of periodic electrical activity, in spite of wave breaks.

The instability in spiral waves occurs because spatial gradients in k (Aliev-Panfilov) or in G (ORd) induce spatial variations in both and : In our simulation domain, the local value of () decreases (increases) from the top to the bottom. In the presence of a single spiral wave (left panel of Fig 4), the domain is paced, in effect, at the frequency ω of the spiral, i.e., with a fixed time period T = 1/ω = APD+ DI, where DI is the diastolic interval (the time between the repolarization of one AP and the initiation of the next AP). Thus, the bottom region, with a long APD, has a short DI and vice versa. The restitution of the conduction velocity θ implies that a small DI leads to a low value of θ and vice versa [40] (see Fig 5). To compensate for this reduction of θ, the spiral wave must reduce its wavelength λ, in the bottom, large-APD (small-DI) region, so that its rotation frequency remains unchanged, as shown in Fig 4 (also see S2 Video), where the shortening of the spiral arms is indicated by the variation of λ along the spiral arm (λ2 > λ1, in the pseudocolor plot of V in the top-left panel t = 1.46 s). Clearly, this shortening is anisotropic, because of the uni-directional variation in k or G; this anisotropy creates functional heterogeneity in wave propagation, causing a local conduction block, which leads in turn to the spiral-wave instability we have discussed above (Fig 4). The phenomenon of conduction block in a medium with a gradient in ionic properties has been extensively investigated in an earlier study [41]; here, it is this local conduction block (caused by the anisotropy of the medium) that leads to the break-up of the spiral arms. It should be noted that the stability of the spiral wave depends on the APD difference between the region, where the spiral is initiated, and the top region, where the APD is maxmium; therefore, its stability depends on the location of the spiral-wave initiation along the vertical direction.

Fig 4

Mechanism of spiral-wave instability.

Pseudocolor plots of V illustrating the development of a spiral-wave instability, with the passage of time t, in the 2D ORd model, with a spatial gradient in G. The left frame shows the shortening of the spiral arm (λ varies along the spiral arm, and λ2 > λ1) indicated just before the spiral wave breaks (see the middle and the right frames).

Fig 5

Cunduction-velocity (CV) restitution curve.

Figure showing a CV restitution curve, generated by using the ORd model. The value of conduction velocity θ initially increases with the increase of diastolic interval DI and saturates at large values of DI.

In the ORd model, we find that gradients in G easily induce instabilities of the spiral for small values of Δg ≡ g − g ≃ 0.5; by contrast, in a medium with gradients in G, the spiral remains stable for values of Δg as large as 4.8 (shown in Fig 6). This implies that the stability of the spiral depends on the magnitude of the gradient in ω that is induced in the medium.

Fig 6

Dependence of spiral-wave stability on Δω.

Pseudocolor plots of V illustrating the development of a spiral-wave instability for (top-row panels) and the sustenance of a stable spiral for (bottom-row panels). The parameter sets are as follows: top row, Δg = 0.5, g = 0.5 and G = G; bottom row, Δg = 0.5, g = 0.2 and G = G.

Scroll-wave instability

In Fig 7 (also see S3 Video), we extend our study to illustrate the onset of scroll-wave instabilities in a 3D, anatomically realistic human-ventricular domain, in the presence of spatial gradients in G. In mammalian hearts, the APD is typically lower in the apical region as compared to that in the basal region [16]. Therefore, we use values of the APD that increase from the apex to the base (and, hence, ω decreases from the apex to base). With g(G) = 6 and Δg = 4, we observe breakup in a scroll wave that is otherwise stable in the absence of this spatial gradient. We note that the mechanism for the onset of such scroll-wave instabilities is the same as in 2D, and it relies on the gradient-induced anisotropic shortening of the scroll wavelength. For control, we also perform a simulation with small Δg = 0.1 that does not show scroll-wave instability (see S4 Video).

Fig 7

Scroll-wave instability.

Pseudocolor plots of V illustrating the occurrence of scroll-wave instability in our anatomically realistic human-ventricular domain, in the presence of an apico-basal gradient in G.

Discussion

We have shown that gradients in parameters that affect the APD of the constituent cells induce spatial gradients in the local value of ω. This gradient in the value of ω leads to an anisotropic reduction in the wavelength of the waves, because of the conduction-velocity restitution property of the tissue, and it paves the way for spiral- and scroll-wave instability in the domain. We would like to point out that this instability is not because of the condition of steep APD restitution curves as reported in ref. [12]. We find that the value of slopes of APD restitution curves for all values of g for all G in Fig 2 are less than one. Therefore, the instability of waves in our study is induced by the anisotropic variation of APD in the medium. This gradient-induced instability is a generic phenomenon because we obtain this instability in the simple Aliev-Panfilov and the detailed ORd model for cardiac tissue. Such an instability should be observable in any excitable medium that has the conduction-velocity-restitution property. We find that the spiral or scroll waves always break up in the low-ω region. This finding is in line with that of the experimental study by Campbell, et al., [15] on neonatal-rat-ventricular cell cultures and a computational study by Xie, et al., [42], who observe spiral-wave break-up in regions with a large APD. We find that the stability of the spiral is determined by the magnitude of the gradient in ω; the larger the magnitude of the gradient in the local value of ω, the more likely is the break up of the spiral or scroll wave. By using the ORd model, we find that ω varies most when we change G (as compared to other ion-channel conductances) and, therefore, spiral waves are most unstable in the presence of a gradient of G. By contrast, we find that ω varies most gradually with G, and hence the spiral wave is most stable in the presence of a gradient in G (as compared to gradients in other conductances).

Earlier studies have investigated the effects of ionic-heterogeneity on spiral-wave dynamics. The existence of regional ionic heterogeneities have been found to initiate spiral waves [43], attract spiral waves to the heterogeneity [44], and destabilize spiral waves [45]. The presence of APD gradients in cardiac tissue has been shown to drive spirals towards large-APD (low-ω) regions [46] or small- APD regions [47], called ‘anomalous drift’, by varying model parameters. We have also observed the drift of spiral waves towards the large-APD region (see Fig 8) in the initial time before the waves break up. A study by Zimik, et al., [35] finds that spatial gradients in ω, induced by gradients in the density of fibroblasts, can precipitate a spiral-wave instability. However, none of these studies provides a clear understanding of the mechanisms underlying the onset of spiral- and scroll-wave instabilities, from a fundamental standpoint. Moreover, none of these studies has carried out a detailed calculation of the pristine effects of each individual major ionic currents, present in a myocyte, on the spiral-wave frequency; nor have they investigated, in a controlled manner, how gradients in ion-channel conductances lead to spiral- or scroll-wave instabilities. Our work makes up for these lacunae and leads to specific predictions that should be tested experimentally.

Fig 8

Spiral-wave drift towards low-ω region.

Figure shows the loci of a spiral tip, which was tracked using the method in ref. [49]. Time is shown in the vertical axis and the X and Y axes are the spatial dimensions of our simulation domain. The gradient in ω is induced by applying a linear gradient in G along the Y axis (indicated by the black arrow), and the spiral drifts along the gradient as indicated by the red arrow. The g and g values are 3 and 1, respectively.

We end our paper by discussing certain limitations of our work. We have shown that large spatial gradients in APD can induce scroll-wave breaks in real hearts via a representative simulation on anatomically realistic heart domain with fiber orientation; however, we have not incorporated other important physiological details of real mammalian hearts, like the intrinsic heterogeneities that exists in them [16, 27–31], and the bidomain nature of the tissue [48]. Moreover, in our study we induce heterogeneity in the medium by applying a spatial gradient that extends throughout the domain, but, heterogeneities in real hearts tend to occur in localized regions. However, our results of spiral-wave break at large-APD region should still hold even if the heterogeneites are localized, as has been shown in [42].

Spiral-wave instability in the Aliev-Panfilov model.

Video of pseudocolor plots of transmembrane potential V showing the formation of spiral-wave instability in a medium with gradient in k: g = 0.5 and g = 1.5. For the video, we use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time.

(AVI)

Click here for additional data file.

Spiral-wave instability in the ORd model.

Video pseudocolor plots of transmembrane potential V showing the formation of spiral-wave instability in a medium with a gradient in G (g = 0.2 and g = 2). For the video, we use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time.

(AVI)

Click here for additional data file.

Video pseudocolor plots of transmembrane potential V showing the formation of scroll-wave instability in an anatomically realistic model for human ventricles. A linear gradient in G is applied along the apico-basal direction: g = 2 in the apex and g = 6 in the base. For the video, we use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time.

(AVI)

Click here for additional data file.

Stable scroll-wave.

Video pseudocolor plots of transmembrane potential V showing a stable scroll-wave for small Δg = 0.1 in an anatomically realistic model for human ventricles. A linear gradient in G is applied along the apico-basal direction: g = 2 in the apex and g = 2.1 in the base. For the video, we use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time.

(AVI)

Click here for additional data file.