Instability of spiral and scroll waves in the presence of a gradient in the fibroblast density: the effects of fibroblast-myocyte coupling

Fibroblast-myocyte coupling can modulate electrical-wave dynamics in cardiac tissue. In diseased hearts, the distribution of fibroblasts is heterogeneous, so there can be gradients in the fibroblast density (henceforth we call this GFD) especially from highly injured regions, like infarcted or ischemic zones, to less-wounded regions of the tissue. Fibrotic hearts are known to be prone to arrhythmias, so it is important to understand the effects of GFD in the formation and sustenance of arrhythmic re- entrant waves, like spiral or scroll waves. Therefore, we investigate the effects of GFD on the stability of spiral and scroll waves of electrical activation in a state-of-the- art mathematical model for cardiac tissue in which we also include fibroblasts. By introducing GFD in controlled ways, we show that spiral and scroll waves can be unstable in the presence of GFDs because of regions with varying spiral or scroll-wave frequency {\omega}, induced by the GFD. We examine the effects of the resting membrane potential of the fibroblast and the number of fibroblasts attached to the myocytes on the stability of these waves. Finally, we show that the presence of GFDs can lead to the formation of spiral waves at high-frequency pacing.


Introduction
The mechanical contractions of the heart muscles are mediated by electrical waves of activation in cardiac tissue. Disturbances in the normal propagation of these electrical waves can be arryhthmogenic because of the excitation of pathological re-entrant waves, such as spiral waves, in two dimensions (2D), and scroll waves, in three dimensions (3D). Spiral waves are linked to cardiac arrhythmias, such as ventricular tachycardia (VT) and ventricular fibrillation (VF). Ventricular tachycardia is known to be associated with the existence of a single spiral (or scroll) wave in the medium [1,2]; life-threatening VF is associated with multiple or broken spiral (or scroll) waves [3,4,5]. Some studies indicate that episodes of VF are preceded by VT [6,7], suggesting the possibility of an initial onset of a spiral wave, which then degenerates to a multiple-spiral-wave turbulent state (VF) [8,9]. Given that VF can prove to be fatal, it is important to understand this transition from VT to VF. Many factors can affect the stability of a spiral wave; hence multiple mechanisms have been proposed for this VT-VF transition [10]. Here, we discuss the effects of inhomogeneities in the fibroblast density distribution on the stability of spiral and scroll waves of electrical activation in a stateof-the-art mathematical model for cardiac tissue, in which we also include fibroblasts.
Fibroblast proliferation (fibrosis) occurs in the heart during myocardial remodelling in, e.g., ischemic hearts or hypertensive hearts [11]; such proliferation is part of the wound-healing process after injuries caused, e.g., by myocardial infarction. When fibroblasts are coupled to a cardiac myocyte, they can change the electrophysiological properties of the myocyte [12,13,14,15]. This modulation of the electrophysiological properties of myocytes can, in turn, alter the dynamics of waves of electrical activation in cardiac tissue. Also, because of the heterocellular coupling between fibroblasts and myocytes, the presence of fibroblasts modulates the conduction properties of cardiac tissue. Therefore, the interposition of fibroblasts between myocytes in cardiac tissue can lead to the fragmentation of the electrical waves [16,17] and even waveblock [17]. Many studies have investigated the effects of fibroblasts on wave dynamics in cardiac tissue [16,17,18,19,20]. Some of these studies model the fibroblasts as inexcitable obstacles [19,20]; others take into account the fibroblast-myocyte coupling and consider either (a) a random distribution of fibroblasts, with an average density that is uniform in space [16,17,18], or (b) localized fibroblast inhomogeneities [18]. However, in real diseased hearts the distribution of fibroblasts may not be uniform, even on average, but, rather, there may be a gradient in fibroblast density (GFD), as has been observed in aged-rabbit hearts [6]. Moreover, in hearts that have been injured, say because of myocardial infarction, the fibroblast density may vary from a high value in the infarcted region to a lower value in the normal region of the heart [21,22], with intermediate values in interfaces between these regions. It is important, therefore, to understand what role such GFD can play in inducing and then, perhaps, destabilizing re-entrant waves, like spiral or scroll waves.
We show that a state-of-the-art mathematical model for cardiac tissue, based on the O'Hara-Rudy model (ORd) for a human ventricular cell [23], provides us with a natural platform for (a) incorporating fibroblast-myocyte interactions and (b) imposing GFD in a controlled way so that we can study, exclusively, its effects on spiral-and scroll-wave dynamics, without other complicating factors that can be present in real cardiac tissue, such as scars, which lead to conduction inhomogeneities [24]. We carry out such a controlled study of the effects of GFD by using the ORd model, for a human ventricular cell [23], with passive fibroblasts, as in the model of MacCannell, et al. [15], which interact with myocytes (see Materials and Methods). We introduce fibroblasts in such a way that we can control GFD systematically. Before we present the details of our study, we give a qualitative overview of our principal results.
We find that GFD can destabilize spiral or scroll waves, i.e., lead to the break up of a spiral or scroll wave into multiple waves. Our in silico studies help us to uncover the principal mechanism of this GFD-induced spiral-or scroll-wave instability. The myocyte-fibroblast coupling changes the action potential duration (APD) of the attached myocytes. We show first how the spiral-wave frequency ω, in a homogeneous domain with randomly distributed fibroblasts, depends inversely on the APD of the myocytes in the medium, a result that is consistent with dimensional analysis. Roughly speaking, GFD induces a gradient in the mean APD, i.e., APD, in the medium; therefore, if a spiral wave forms, with its core in a high-APD region, its rotation frequency near the core is low. Such a low-ω region cannot support wave trains that come from a high-ω region, say because of a spiral wave there. We find that the greater the variation of APD with GFD the more readily does spiral-or scroll-wave instability set in. We build on this qualitatively appealing argument by carrying out detailed numerical investigations of GFD-induced spiral-and scroll-wave break up, in two-and three-dimensional simulation domains, respectively, and in three models, the first with a uniform distribution of fibroblasts, and two others in which the mean density of fibroblasts changes either (a) linearly along one spatial direction or (b) as a step function. Furthermore, we investigate the effects of other fibroblast parameters, like the resting membrane potential E f and the number N f of fibroblasts attached to myocytes, on the stability of such spiral and scroll waves. For a given GFD, we find that the larger the values of E f and N f , the more readily does spiral-or scroll-wave instability set in. Finally, we show how the presence of GFD can lead to the formation of re-entrant waves, like spiral waves, when we pace an edge of our simulation domain at a high frequency.
The remaining part of this paper is organized as follows. The Materials and Methods Section describes the models we use and the numerical methods we employ to study them. The Section with the title Results contains our results, from single-cell and tissue-level simulations. Finally, the Section called Discussions is devoted to a discussion of our results in the context of earlier numerical and experimental studies.

Materials and Methods
We use the state-of-the-art O'Hara-Rudy model (ORd) for a human ventricular cell [23] for our myocyte cell. In the ORd model, we implement the modifications suggested in Ref. [23], where the fast sodium current (I Na ), of the original ORd model, has been replaced with that of the model due to Ten Tusscher and Panfilov [25]. The fibroblasts are modelled as passive cells, for which we use the model given by MacCannell, et al. [15], but, instead of using a constant membrane conductance G f , we use a G f that has a nonlinear dependence on its membrane voltage V f , as in Refs. [16,26]. The value of G f is set to 1 nS if V f is below -20 mV, and 2 nS for V f above -20 mV. The gap-junctional conductance between a myocyte and a fibroblast is set to 8 nS.
In our two-dimensional (2D) simulations the fibroblasts are attached atop the myocytes; thus, our 2D simulation domain is a bilayer, as in previous studies [18,27]; similarly, we place fibroblasts in the interstitial region between myocytes in our threedimensional (3D) simulations as in Ref. [28]. The number of fibroblasts N f attached to a myocyte in a fibroblast-myocyte composite is 2, unless otherwise mentioned in the text. For a given percentage P f of fibroblasts in the medium, the fibroblasts are attached to the myocytes in a tissue with a probability P f /100. We carry out three types of simulations for three different models of fibroblast distribution: Model-I, Model-II, and Model-III. In Model I, the value of P f (x, y) is constant throughout the domain, as in equation (1) (see figure 2 (a)). Thus, the distribution of fibroblast is uniform and isotropic. In Model II, the value of P f (x, y) varies linearly along the vertical direction (y axis) of the domain (see figure 2 (b)) given by equation (2). In the last Model III, the value of P f (x, y) changes discontinuously (see figure 2 (c)) given by equation (3).
where p f is a constant between 0 and 100. Model where L is the size of the square domain, and p f1 and p f2 are constants between 0 and 100. Model where Θ(y) is the Heaviside step function, and y 0 is a particular value of the distance along the vertical direction; we use y 0 = 3 4 L in our study. The membrane potential of a myocyte V m is governed by the ordinary differential where C m is the myocyte capacitance, which has a value of 185 pF; I ion is the sum of all the ionic currents of the myocyte, I gap is the gap-juctional current between the fibroblast and the myocyte, and N f is the number of fibroblasts attached to the myocyte. We give I ion and I gap below: we list the ionic currents of the myocyte in Table 1.
The membrane potential of a fibroblast is governed by the equation where C f =6.3 pF is the membrane capacitance of the fibroblast, and I f is the fibroblast current, here E f is the resting membrane potential of the fibroblast, and its range of values in our study is from 0 to -50 mV, as observed experimentally in Refs. [26,29,30,31]. inward rectifier K + current I NaCa Na + /Ca 2+ exchange current I NaK Na + /K + ATPase current I Nab Na + background current I Cab Ca 2+ background current I pCa sarcolemmal Ca 2+ pump current I Kb K + background current I CaNa Na + current through the L-type Ca 2+ channel I CaK K + current through the L-type Ca 2+ channel A list of the ionic currents in the ORd model (symbols as in Ref. [23]).
The spatiotemporal evolution of the membrane potential (V m ) of the myocytes in tissue is governed by a reaction-diffusion equation, which is the following partialdifferential equation (PDE): where D is the diffusion constant between the myocytes.

Numerical Methods
We solve the ODEs (4) and (7) for V m and V f , respectively, and also the ODEs for the gating variables of the ionic currents of the myocyte by using a forward-Euler method. For solving the PDE (9), we use the forward-Euler method for time marching with a five-point stencil for the Laplacian in 2D and a 7-point stencil in 3D. We set D= 0.0012 cm 2 /ms. The temporal and spatial resolutions are set to be δx= 0.02 cm and δt=0.02 ms, respectively. The conduction velocity of a plane wave in the tissue, with the above set of parameters, is 65 cm/s. In our two-dimensional (2D) tissue simulations, we use a domain size of 960 × 960 grid points, which translates into a physical size of 19.2 × 19.2cm 2 . We initiate the spiral wave by using the conventional S1-S2 cross-field protocol, where we first apply an S1 plane wave and allow its wave back to cross some part of the domain (see figure 1, time= 300 ms) and then we apply the S2 stimulus perpendicular to the S1 wave (figure 1, time= 320 ms), as shown in figure 1. The strength and duration of both the S1 and S2 stimuli are -150 µA/µF and 3 ms, respectively. All our tissue simulations are carried out for a duration of 10 seconds. . S1-S2 cross-field protocol. Pseudocolor plot illustrating the S1-S2 cross-field protocol that we use to initiate a spiral wave in the domain. The colorbar indicates the membrane voltage V m in millivolts.

Results
We first study the effects of the gap-junctional coupling G gap of fibroblasts on the electrophysiological properties of myocytes. We then investigate the effects of G gap on the spiral-wave frequency in a 2D simulation. Next we show how GFD leads to spiral-wave instability. Thereafter, we show the effects of GFD on scroll-wave stability in a 3D simulation domain. Finally, we discuss the role of GFD on the initiation, via pacing, of re-entrant spiral waves.

Effects of fibroblast gap-junctional coupling, on a myocyte, and the spiral-wave frequency
The gap-junctional coupling of fibroblasts with a myocyte changes the electrophysiological properties of the myocyte [12,13,14,15]. For instance, it can modulate the action potential duration (APD) of a myocyte. The APD of the myocyte may increase or decrease depending on the electrophysiological properties of the fibroblasts, like their resting-membrane potential E f . Figure 3 shows the action potentials of a myocyte attached to N f = 2 fibroblasts with different values of E f . The APDs of the myocyte attached to fibroblasts, with E f = -15 mV (blue curve) and E f = -25 mV (black curve) are larger than that of an isolated myocyte (red curve), whereas the APD of the myocyte attached to fibroblasts with E f = -50 mV (magenta curve) is lower than that of an isolated myocyte. This fibroblast-induced modulation of the APDs of the myocytes affects the properties of electrical waves, like the spiral-wave frequency, at the tissue level. This spiral frequency ω depends on the APDs of the constituent myocytes as follows. Consider a stable spiral that does not meander; dimensional analysis yields where θ is the conduction velocity and λ is the wavelength. Furthermore, λ ≃ θ × APD, and, therefore, In a domain with a uniform and isotropic fibroblast distribution with P f (x, y)=p f (Model-I, equation 1 and figure 2 (a)), the average APD (APD) of the myocytes depends on p f , and, therefore, the frequency ω of a spiral wave depends on p f . In figure 5 we plot ω versus p f for E f = -50mv (red curve), and E f = -25 mV (blue curve), and E f = -15 mV (black curve). We see from the plot for E f = -15 mV that ω decreases as p f increases; this decrease is slower for E f = -25 mV; and, by the time E f = -50 mV, ω increases with p f . The decrease of ω with p f for E f = -25 mV and -15 mV is because fibroblasts, with these values of E f , increase the APD of an attached myocyte (see figure 3); and, therefore, in a tissue with a uniform and isotropic fibroblast distribution, the increase of p f increases APD of the constituent myocytes, and, thereby, decreases ω (see equation (11)). However, for E f = -50 mV, the APD of the myocyte decreases (see figure 3), so ω increases with p f . For the plots in figure 5 we obtain the spiral frequency ω from the averaged power spectra of time-series of V m , which we record from four representative points of the domain (white squares in figure 4 (a)). We show one representative example for p f = 40% in figure 4. Figure 4(a) shows the spiral wave in the domain; and figure 4(b) shows the power spectrum of V m averaged over those from the points (white squares in figure 4(a)) of the domain. The spiral-wave frequency ω is the value of the dominant peak in the spectrum, which is 4.4 Hz in this case. No fibroblast E f = −50mV , and E f = -50 mV (magenta curve). The APDs of the myocyte attached to fibroblasts of E f = -15 mV and E f = -25 mV are larger than that of the isolated myocyte, whereas the APD of the myocyte attached to fibroblasts of E f =-50 mV is lower than that of the isolated myocyte.

Instability of spiral waves
The spiral-wave frequency ω, in a domain with fibroblasts, depends on the fibroblast density p f . Therefore, in a domain with an inhomogeneous distribution of fibroblasts, the spiral-wave frequency (the frequency at which the spiral tip rotates) varies with space. To illustrate this, we take the Model-II fibroblast distribution (see equation 2), where P f (x, y) varies along the y axis. This variation of P f (x, y) in the y direction induces a y dependence in APD and, therefore, in ω. We set p f1 = 10%, p f2 = 80%, and E f = -25 mV. Figure 6 show the cases when a spiral wave is initiated in three different regions of the domain: (a) top, (b) middle, and (c) bottom. The frequencies ω of the spiral tip in these top, middle, and bottom regions are 4.2 Hz, 4.3 Hz, and 4.4 Hz, respectively. We measure these frequencies from the time series of V m , which we record from the four representative points, shown in figure 4 (a). The spatial dependence of the local value of ω in the domain can lead to spiral-wave instability if the variation in ω is sufficiently large. We show this for the Model-II fibroblast distribution, with p f1 = 10%, p f2 = 100%, and E f = -25 mV. If we initiate a spiral somewhere in the middle of our simulation domain, it breaks up as shown in figure 7 in the upper region of the domain, where ω is low (see also the Supplementary Video S1). This can be understood qualitatively as follows: The upper, low-ω region cannot support the high-frequency wave-trains that are emitted by the spiral tip, which is rotating in the middle region, with a higher value of ω than in the upper part of the domain; the inability of the upper region to support high-frequency waves leads to anisotropic conduction blocks because of the anisotropic fibroblast distribution and, hence, wave breaks.
This breaking of waves in the low-ω region (i.e., large APD region) has also been observed experimentally in monolayers of neonatal-rat ventricular myocytes by Campbell, et al. [32]. In these experiments, the gradient in APD has been induced by varying the I Kr ion-channel density.
We have shown that spiral-wave instability stems from the spatial gradient in ω that is induced by the GFD. We might expect, therefore, that steep changes in ω might lead to an enhancement in such instability. To show this, we study the following two cases: case (A), with the Model-II fibroblast distribution, p f1 = 10%, and p f2 = 75%; and case (B), with the Model-III fibroblast distribution, given in equation 3 and illustrated in figure 2 (c), with p f1 = 10% and p f2 = 65%. The value of E f for both these cases is -25 mV. We show in figure 8   higher in case (A), with the Model-II distribution, than in case (B), with the Model-II distribution (see also the Supplementary Video S2). Therefore, a comparison of our results for cases (A) and (B) shows that a high local slope of dω dy , which has a singularity at y= 3 4 L in Model-III, enhances spiral-wave instability.

Dependence of spiral-wave instability on fibroblast parameters
The APD of a myocyte, in a fibroblast-myocyte composite, depends on E f , and, therefore, the stability of a spiral wave in a domain with GFD also depends on E f . To illustrate this, we first recall the variation of ω with p f for three values of E f in figure 5, namely, E f = -15 mV (black curve), E f = -25 mV (blue curve), and E f = -50 mV (red curve). Note that the change in ω with p f is more significant for E f = -15 mV than for E f = -25 mV. Therefore, in a domain with a given GFD, the spiral breaks up with  The variation of ω with p f also depends on the number N f of fibroblasts that are attached to the myocytes. Figure 11 shows the variation of ω with p f for different values of N f , for a constant value of E f = -15 mV. The black, blue, and red curves are for N f = 1, 2, and 3, respectively. These plots show that the larger the value of N f the larger is the change of ω with p f . Hence, in a domain with a given GFD and E f , the spiral-wave instability increases with N f , because the variation of ω with p f increases with N f . To illustrate this, figure 12 (top panel) shows the breaking of a spiral wave for a GFD with the Model-II distribution, p f1 = 10%, p f2 = 60%, E f = -15 mV, and N f = 3. The same GFD and E f does not lead to spiral break-up for N f = 2 ( figure 12, bottom panel). (See also the Supplementary Video S5.)

Three-dimensional simulation domains
To check if our results also hold in 3D simulation domains, we perform a few representative 3D simulations, where the thickness (x dimension) is 2 mm and the linear dimensions in the y and z directions are both 19.2 cm. We initiate a scroll wave by using the same S1-S2 cross-field protocol that we use in our 2D simulations for spiral-wave initiation. We find that the results we obtain in 3D, for the dependence of scroll-wave stability on E f and N f , are qualitatively similar to those we have obtained above for spiral-wave stability in 2D. For example, we show the dependence of scroll-wave stability  on N f . We take a 3D domaim with GFD given by the Model-II fibroblast distribution: P f (x, y, z) varies from p f1 = 10% to p f2 = 60% along the z axis, but is constant in the x-y plane for a fixed value of z; we use E f =-15 mV here. Now we study the stability of a scroll wave for N f = 2 and N f = 3. Figure 13(a) shows a scroll wave breaking up for N f = 3; however, the scroll does not break up for N f = 2, as shown in figure 13(b) (see also the Supplementary Video S6). Therefore, with GFD, the higher the value of N f the greater is the instability of scroll waves (as in our 2D simulations).

Pacing-induced re-entry
So far we have discussed the stability of spiral and scroll waves in the presence of GFD, where we initiate the spiral or scroll waves by using the S1-S2 cross-field protocol. We now show GFD itself can lead to re-entry (i.e., the formation of spiral or scroll waves) if we pace the system. To illustrate this, we pace our 2D simulation at a high frequency, with a pacing cycle length PCL= 250 ms, in the presence of GFD (the Model-II fibroblast distribution, p f1 = 10%, p f2 = 100%, and E f = -15 mV). The pacing stimulus is applied at the bottom edge of the domain. We apply 20 pulses with this PCL, and we observe that a spiral wave is induced by such high-frequency pacing, as shown in figure 14 (see also the Supplementary Video S7). Thus, GFD not only provides a substrate for spiral-or scroll-wave instabilities in cardiac tissue, but it can also be a source of such re-entrant  waves if the tissue is paced at sufficiently high frequency. Such pacing-induced spiral waves occur with high-frequency pacing (low PCL) but not with low-frequency (large PCL) pacing. In the stability diagram of figure 15, we show the region in which we observe such pacing-induced re-entry in the E f -PCL plane in the presence of the GFD. The black-colored region indicates the region in which we obtain spiral waves, and the red-colored region is the region in which we do not see spiral waves. Thus, for a given E f , re-entry occurs below a certain PCL; and, as E f decreases, the threshold of PCL, below which re-entry occurs, decreases. Re-entry occurs with high-frequency and not with low-frequency pacing for the following reasons. The wavebacks of the travelling waves repolarize (come back to the resting state) slowly in the low-ω regions, which are in the top part of our simulation domain in figure 14, because APD is large in these regions. When we use low-frequency pacing, the wave-fronts do not meet the wavebacks of the preceding waves. By contrast, when we use high-frequency pacing, the wavefronts of the succeeding waves interact with the wavebacks of the preceding waves, and, because of the anisotropic repolarization that arises from the GFD, the succeeding wavefronts become corrugated, which leads to re-entry and the formation of spiral waves.

Discussion
We have shown by detailed numerical simulations how gradients in the fibroblast density (GFD) affect spiral-or scroll-wave dynamics in the state-of-the-art ORd mathematical model for cardiac tissue. Our work shows that such gradients induce spatial variations in the local spiral-or scroll-wave frequency ω. These variations play a crucial role in the stability of these spiral or scroll waves in the presence of GFD. In particular, we find that spiral or scroll waves break in regions where the local ω is low. Furthermore, for a particular GFD, the stability of these waves depend on E f and N f , insofar as they affect the variation of ω with p f (see figures 5 and 11). Last, but not least, high-frequency pacing at one end of our domain can lead to the formation of spiral waves if GFD is present.
Some earlier studies have investigated the effects of fibroblasts on spiral-and scroll-wave dynamics. Numerical simulations of mathematical models of cardiac tissue, with a distribution of fibroblasts, have shown single or multiple spiral or scroll waves [17], depending on the fibroblast density. A detailed study of the effects Stability diagram in the E f -PCL plane. Figure showing the reentry (black) and no re-entry (red) regions in the E f -PCL plane (see text). This figure shows that re-entry occurs in the low-PCL region; and, as E f decreases, the threshold PCL, below which re-entry occurs, also decreases.
of homogeneous or localized distributions of fibroblasts on spiral-wave dynamics has been carried out by Nayak, et al. [18]. Fibroblasts, randomly distributed, have been shown to convert complex wave patterns, like mutiple-spiral states, to simpler patterns, like a single rotating spiral state [28]. In a mathematical model of heartfailure tissue, Gomez, et al., have found that the presence of fibroblasts enhances re-entrant phenomena [33]. Also, the fibroblast-myocyte coupling has been found to be arrythmogenic, because of its ability to induce pathological excitations, like alternans [12] and early afterdepolarizations [6]; the latter can, in turn, lead to reentries.
Although these studies, and many others [19,20,21], have investigated the causes of re-entrant waves and their sustenance in cardiac tissue with fibroblasts, none of them has studied, systematically, spiral-and scroll-wave instability because of the fibroblastmyocyte coupling in a domain with GFD. Our study, in which we control GFD, has allowed us to study clearly the arrhythmogenic potential of such GFD. It is important to carry out such an investigation because the gradients in density distribution of fibroblasts are known to occur in real hearts [6]. We hope our results will lead to detailed studies of GFD-induced spiral-and scroll-wave instability at least in in vitro experiments on cell-cultures. At the simplest level, we suggest fibroblast analogs of the experiments of Campbell, et al. [32], in which the gradient in APD has been induced by varying the I Kr ion-channel density.
We end our discussion with some limitations of our study. Our tissue simulations do not take into account anisotropy because of the orientation of muscle fibers. Such tissue anisotropy may exacerbate the instability of spiral waves in the presence of a GFD. We have used a monodomain representation of cardiac tissue; our study needs to be extended to other tissue models, such as those that use bidomain representations [34]. We use a passive model of the fibroblasts; however, there is some evidence [15,35,36], that fibroblasts can behave as active cells. Nonetheless, our qualitative results, about the arrhythmogenic effects of GFD, and its root cause, should not depend on such details.
in the bottom part of the domain because it is the low-ω region. We use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time. (MPEG) Video S5. Video showing the break-up of a spiral in the presence of GFD with the Model-II distribution, p f1 =10%, and p f2 =60%, for E f = -15 mV and N f = 3 (left panel). The same GFD and E f does not lead to break-up for N f = 2 (right panel). We use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time. (MPEG) Video S6. Video (bottom panel) showing the break-up of a scroll wave in a 3D domain with GFD (Model-II distribution, with P f (x, y, z) varying from p f1 = 10% to p f2 = 60% along the z-axis and constant in the x-y plane for a fixed value of z), where the value of E f is set to -15 mV and N f =3. The same GFD and E f does not lead to break up for N f =2 (top panel). We use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time. (MPEG) Video S7. Video showing the pacing-induced occurrence of a spiral wave, in a 2D simulation domain with GFD (the Model-II distribution, p f1 = 10% and p f2 = 100% GFD, and E f = -25 mV) paced with PCL= 250 ms. We use 10 frames per second with each frame separated from the succeeding frame by 20ms in real time. (MPEG)