Spiral-Wave Turbulence and Its Control in the Presence of Inhomogeneities in Four Mathematical Models of Cardiac Tissue

Regular electrical activation waves in cardiac tissue lead to the rhythmic contraction and expansion of the heart that ensures blood supply to the whole body. Irregularities in the propagation of these activation waves can result in cardiac arrhythmias, like ventricular tachycardia (VT) and ventricular fibrillation (VF), which are major causes of death in the industrialised world. Indeed there is growing consensus that spiral or scroll waves of electrical activation in cardiac tissue are associated with VT, whereas, when these waves break to yield spiral- or scroll-wave turbulence, VT develops into life-threatening VF: in the absence of medical intervention, this makes the heart incapable of pumping blood and a patient dies in roughly two-and-a-half minutes after the initiation of VF. Thus studies of spiral- and scroll-wave dynamics in cardiac tissue pose important challenges for in vivo and in vitro experimental studies and for in silico numerical studies of mathematical models for cardiac tissue. A major goal here is to develop low-amplitude defibrillation schemes for the elimination of VT and VF, especially in the presence of inhomogeneities that occur commonly in cardiac tissue. We present a detailed and systematic study of spiral- and scroll-wave turbulence and spatiotemporal chaos in four mathematical models for cardiac tissue, namely, the Panfilov, Luo-Rudy phase 1 (LRI), reduced Priebe-Beuckelmann (RPB) models, and the model of ten Tusscher, Noble, Noble, and Panfilov (TNNP). In particular, we use extensive numerical simulations to elucidate the interaction of spiral and scroll waves in these models with conduction and ionic inhomogeneities; we also examine the suppression of spiral- and scroll-wave turbulence by low-amplitude control pulses. Our central qualitative result is that, in all these models, the dynamics of such spiral waves depends very sensitively on such inhomogeneities. We also study two types of control schemes that have been suggested for the control of spiral turbulence, via low amplitude current pulses, in such mathematical models for cardiac tissue; our investigations here are designed to examine the efficacy of such control schemes in the presence of inhomogeneities. We find that a local pulsing scheme does not suppress spiral turbulence in the presence of inhomogeneities; but a scheme that uses control pulses on a spatially extended mesh is more successful in the elimination of spiral turbulence. We discuss the theoretical and experimental implications of our study that have a direct bearing on defibrillation, the control of life-threatening cardiac arrhythmias such as ventricular fibrillation.


Introduction
Cardiac arrhythmias like ventricular tachycardia (VT) and ventricular fibrillation (VF) are a major cause of death in industrialised countries. Experimental studies over the past decade or so have suggested that VT is associated with an unbroken spiral wave of electrical activation on cardiac tissue but VF with broken spiral waves [1][2][3]. There is growing consensus [4,5] that the analogs of VT and VF in mathematical models for cardiac tissue are, respectively, (a) a single rotating spiral wave in two dimensions or a scroll wave in three dimensions and (b) spiral-wave or scrollwave turbulence [6][7][8]. It is imperative, therefore, to study spiraland scroll-wave turbulence systematically in such mathematical models. There are several such studies [9][10][11] but none, to the best of our knowledge that compares several models with a view to elucidating low-amplitude defibrillation schemes, which are designed to eliminate spiral-wave turbulence, especially in the presence of inhomogeneities, such as conduction inhomogeneities in cardiac tissue. We address this question here by considering four models of cardiac tissue that are, in order of increasing complexity, (a) the Panfilov model [12], (b) the Luo-Rudy Phase I (LRI) model [13], (c) the reduced Priebe-Beuckelmann (RPB) model [14], and (d) the TNNP model of ten Tusscher, et al. [15]. The Panfilov model [12] is of the Fitzhugh-Nagumo [16,17] type with two fields, namely, the transmembrane potential V and the recovery variable g, which depend on space and time; it is much simpler than models that account for ion channels in the membrane; nevertheless, it yields spiral waves and their break up in a manner that is qualitatively similar to such pattern formation in more realistic ionic models; given its simplicity, the Panfilov model is a good testing ground for numerical and semi-analytical studies. Realistic ionic models, such as the LRI, RPB, and TNNP models, are based on the Hodgkin-Huxley [18] formalism for ion channels. The LRI model [13] uses data from measurements on guinea-pig myocardial cells and accounts for ion channels and voltage dependent ion-channel gating variables. The RPB model [14] improves on the LRI model by incorporating data obtained from human ventricular muscles; furthermore, it includes an ion pump and an ion exchanger. The TNNP model is based on recent experimental data from human ventricular cells; it includes an ion pump, an ion exchanger, and more details of the calcium-ion dynamics [15] than the LRI and RPB models.
Our goal is to carry out a comparison of spiral-wave dynamics in these four models of cardiac tissue especially in the presence of different types of inhomogeneities, such as conduction and ionic inhomogeneities; a comparison of such spiral waves in the Priebe Beuckelmann (PB), RPB, TNNP, models and the model of Iyer et al. [19] has been presented in Ref. [20] but without inhomogeneities. We also study schemes for eliminating spiral-wave chaos from the simulation domain in Panfilov, LRI, RPB, and TNNP models in both homogeneous and inhomogeneous cases. The motivation for undertaking this study is that cardiac tissue contains both conduction and ionic heterogeneities. These can be caused, inter alia, by (a) a myocardial infarction that leads to ischemia [21], the subsequent damage or death of the affected cardiac cells, and, in the latter case, the formation of scar tissue that is nonconducting, (b) chronic heart failure, (c) genetic disorders, or (d) the presence of major blood vessels.

Review of Previous Work
Conduction inhomogeneities in cardiac tissue can affect spiral waves in several ways. Experimental studies [22][23][24] have found that such inhomogeneities can anchor a spiral wave or, in some cases, can even eliminate it completely [2]. Studies of the dependence of such anchoring on the size of the obstacle [22][23][24] reveal that the larger the obstacle the more likely is the anchoring; however, even if the obstacle is large, the wave might not attach to it; furthermore, small obstacles can anchor spiral waves, albeit infrequently [24]. Such behaviors have also been seen in numerical simulations of spiral-wave turbulence in models for cardiac tissue: In particular, Xie et al. [25], have studied the dynamics of spiral waves in the LRI model in a two-dimensional (2D) circular domain with a circular hole in the middle: The parameters and initial condition are so chosen that, without the hole, spiral waves break up in the simulation domain. By shrinking the radius of the hole, the system is changed continuously from a 1D ring to 2D tissue with an obstacle, and, finally, to homogeneous 2D tissue [the hole radius is changed from that of the simulation domain (.9.2 cm) to zero]. When the radius of the hole is very large, the system is effectively a 1D ring; the wave just goes around this ring. As the radius of the hole is decreased, the wave appears as a spiral anchored on the hole but rotating around it periodically, if the hole is large. As the hole radius is decreased a transition occurs first to a quasiperiodically rotating spiral wave and, eventually, to spiral-wave break up and spatiotemporal chaos [25] with the spirals not attached to the hole.
ten Tusscher et al. [26] have studied the Panfilov model with nonexcitable cells distributed randomly in it. In particular, they investigate spiral-wave dynamics as a function of the percentage of the simulation domain covered by such nonexcitable cells and find that, when this percentage is high, spiral-wave break up can be suppressed.
A detailed numerical and analytical study of the interaction of excitation waves with a piecewise linear obstacle has been carried out in Ref. [27]. This study finds that, if the excitability of the medium is high, the wave moves around the obstacle boundary, rejoins itself, and then proceeds as if it had not encountered any obstacles in its path. However, if the excitability is low, the two ends of this wavefront are unable to join, so two free ends survive, curl up, and then develop into two spiral waves, which can in turn break up again. In addition, apart from the excitability of the medium and the local curvature of the wave front, the shape of the obstacle also affects the attachment of spiral waves to it. We have carried out a detailed numerical study [28] of spiral-wave dynamics in the presence of conduction inhomogeneities in the Panfilov and LRI models; our study has shown that the dynamics of spiral waves depends very sensitively on the position of a conduction inhomogeneity.

Summary of Our Results
Ionic inhomogeneities, formed say by local modifications of the maximal conductance of calcium ion channels, affect the action potential of a cardiac cell; in particular, the action potential duration (APD) and other time scales, such as the extent of the plateau region and the refractory period [29], are modified by these inhomogeneities and affect spiral wave dynamics in turn. For example, the stability of a spiral wave, in homogeneous, twodimensional cardiac tissue depends on the maximal amplitude of the slow inward calcium current (governed by the conductance G si ) as illustrated by the numerical study of Qu et al. [30] for the LRI model: As they increased G si they first observed a rigidly rotating spiral wave, then one in which the spiral tip meandered quasiperiodically, and eventually chaotically; finally they obtained spiral turbulence with broken spiral waves. Furthermore, the numerical studies of Refs. [28,31] have found that ionic heterogeneities can play an important role in the initiation and break up of spiral waves; and Ref. [28] has presented preliminary studies of the Panfilov-model analog of ionic inhomogeneities.
We consider spiral-wave dynamics in an otherwise homogeneous medium with a square region in which the conduction or ionic parameters are different from their values in the rest of the simulation domain. We find that such inhomogeneities can have dramatic effects on spiral wave dynamics. We have reported earlier that conduction inhomogeneities can act as anchoring sites for spirals, or lead to the complete elimination of spiral waves, or have no effect on spiral-wave break up; which one of these results is obtained depends on the size and position of the conduction inhomogeneity [28]. In this paper we extend our work to ionic inhomogeneities. We find that such inhomogeneities can also result in the elimination of spiral waves; this depends on the position of the inhomogeneity. Here too we find anchored spirals, but with richer dynamics than with conduction inhomogeneities; e.g., we find states with rotating spiral waves that show period-4 and period-5 cycles and also states that show a coexistence of a periodically rotating spiral-wave and chaotic patterns with broken spiral waves, in the region of the ionic inhomogeneity. Lastly we investigate the efficacy of two low-amplitude schemes [32,33] that have been suggested for the control of spiral-wave turbulence in mathematical models for cardiac tissue. In particular, we carry out detailed simulations of such control schemes in the presence of conduction inhomogeneities; our study shows that the elimination of spiral-wave turbulence is considerably more complicated if inhomogeneities are present than if they are not. This paper is organised as follows: In the Section on ''Methods'' we present the models and numerical methods that we use in our study. In the Section on ''Results'' we present our results on studies of spiral-wave dynamics in the presence of conduction and ionic inhomogeneities; we then give an analysis of two different low-amplitude control schemes for the elimination of spiral-wave turbulence in models for cardiac tissue; finally we examine the efficacies of these control schemes in the presence of conduction inhomogeneities. The concluding Section ''Discussion'' contains a summary of the significance of our results. The Supplementary Material S1 contains equations for the ionic models we study and additional figures that give more details about our results.

Methods
The Panfilov model [12,34] comprises two coupled partial differential equations (PDEs) for the transmembrane potential V (denoted by e in Refs. [12,34]) and the slow, recovery variable g, via which this model accounts, approximately, for the effects of the different ion channels in cardiac tissue; both V and g depend on the spatial coordinate x and time t: The fast initiation of the action potential of this model is encoded in the following piecewise linear form of f(V): The physically appropriate parameters are [12,34] V 1 = 0.0026, V 2 = 0.837, C 1 = 20, C 2 = 3, C 3 = 15, a = 0.06 and k = 3. The time scales of the recovery variable are determined by the function e(V,g) that is e 1 for V,V 2 , e 2 for V.V 2 , and e 3 for V,V 1 and g,g 1 ; g 1 = 1.8, e 1 = 0.01, e 2 = 1.0, and e 3 = 0.3; here we deal with dimensionless quantities. To obtain dimensioned quantities [12,34] we define dimensioned time to be 5 ms times dimensionless time, 1 spatial unit to be 1 mm, and the dimensioned value of the conductivity constant D to be 2 cm 2 /s. In some of our studies we vary e 1 to mimic the effects of ionic inhomogeneties. Voltages in this model are scaled such that the resting potential is zero.
Even though the Panfilov model is very simple when compared to the LRI, RPB, and TNNP models, which retain details of ion channels, it captures several essential properties of the spatiotemporal evolution of V [12,[34][35][36]. In particular, the action potential of the Panfilov model contains both the absolute and relative refractory periods seen in more advanced models. Furthermore, the appearance, propagation, and break up of spiral-wave patterns and the ways in which they can be controlled are similar in these models as we will illustrate here.
The LRI model uses the Hodgkin-Huxley formalism for ion channels in a given cell. It accounts for 6 ionic currents (e.g., Na + , K + , and Ca 2+ ) that flow through voltage-gated ion channels; 9 gate variables regulate the transport of ions across the cell membrane [13]; in the quiescent state the concentration differences of the ions, inside and outside the cell, is such that a potential difference .284 mV is induced across this membrane. Electrical stimuli, which raise the potential across the cell membrane above .260 mV, change the conductances of the ion channels and thus yield an action potential with a typical duration of .200 ms. After the initiation of the action potential, there is a refractory period during which a stimulus of the same strength cannot lead to further excitation of that cell. This excitation moves from one cell to another in the LRI model by virtue of a diffusive coupling; thus the transmembrane potential V obeys a reaction-diffusion-type PDE that includes ionic currents (see Section 2 of the Supplementary Material S1 for details); the time evolution and V dependence of these currents are given by 7 coupled ordinary differential equations (ODEs) [13,28].
We also study the Reduced-Priebe-Beuckelmann (RPB) model [14] that is a simplified version of the Priebe-Beuckelmann (PB) model [37]. The original PB model is a second-generation model based on the phase-2 Luo-Rudy model [38] of a guinea-pig ventricular myocyte with currents scaled to fit human-cell data. In the PB model ion concentrations in a cell can vary in time; and it accounts for some ion pumps. Such second-generation models reproduce ionic currents and concentrations in a single cell during electrical activity; however, the large number of variables in these models pose a significant computational challenge, especially in the simulation of arrhythmias in three-dimensional (3D) and even in two-dimensional (2D) domains. This has motivated the development of the reduced PB (RPB) model, in which the variables are reduced from 15 to 6 by a reformulation of some currents and by approximating all intracellular ionic concentrations by suitable constants. Nevertheless, the RPB model retains important properties of human ventricular tissue such as restitution properties, the shape of the action potential (AP), and the change of this shape as a function of major ionic currents. As in the LRI model, cells in the RPB model have a diffusive coupling with each other; the equations for the RPB model are given in Section 2 of the Supplementary Material S1. In particular, the equilibrium voltage across the cardiac cell membrane is 291 mV in the RPB model.
The most realistic model we study is the one introduced recently [15] by ten Tusscher, Noble, Noble, and Panfilov (TNNP). It is based on experimental data obtained from human ventricular cells. The TNNP model allows for variations of intracellular ion concentration, as in other second-generation models, contains 12 ionic currents, 12 gating variables, one ion pump, and an ion exchanger. All major ionic currents are included in the TNNP model, e.g., the fast inward Na + current I Na , the L-type Ca 2+ current I CaL , the transient outward potassium current I to , the slow, potassium, delayed, rectifier current I Ks , the rapid, potassium, delayed, rectifier current I Kr , and the inward, rectifier K + current I K1 . These and other currents and the details of the dynamics of calcium ions are given in Section 2 of the Supplementary Material S1. As in the LRI and RPB models, cells in the TNNP model are coupled diffusively.
Since this model has many variables, numerical simulations of spiral-wave dynamics in it are considerably harder than in the simpler LRI and RPB models. In both the RPB and TNNP models, we can study human epicardial, endocardial, and M cells by a suitable choice of parameters; for the RPB model we use the parameters for M cells and for the TNNP model we use parameters for epicardial cells (the equilibrium voltage across the cardiac cell in the latter is 286.2 mV).
To integrate the Panfilov model PDEs in d spatial dimensions we use the forward-Euler method in time t, with a time step dt = 0.022, and a finite-difference method in space, with step size dx = 0.5, and five-point and seven-point stencils, respectively, for the Laplacian in d = 2 and d = 3 for spatial grids on square or simple-cubic simulation domains with side L mm, i.e., (2L) d grid points. We use a similar forward-Euler method for the LRI PDEs, with dt = 0.01 ms, and a finite-difference method in space, with dx = 0.0225 cm, and a square simulation domain with side L = 90 mm, i.e., 4006400 grid points. We have checked in representative simulations for the LRI model that a Crank Nicholson scheme yields results in agreement with the numerical scheme described above. The simulation schemes that we use for the RPB and TNNP models are similar to the one we use for the LRI model; for the RPB model we use dt = 0.01 ms, dx = 0.0225 cm, and a 5126512 square simulation domain; for the TNNP case we use dt = 0.02 ms, dx = 0.0225 cm, and a 6006600 square simulation domain (i.e., L = 135 mm).
For all the models that we study, we use Neumann (no-flux) boundary conditions. The initial conditions we use will be specified below. For numerical efficiency, these simulations have been carried out on parallel computers with MPI codes that we have developed for all these models.
As suggested in Ref. [39], it is useful to test the accuracy of the numerical scheme used by varying both the time and space steps of integration. We illustrate this for the TNNP model by measuring the conduction velocity (CV) of a plane wave, which is injected into the medium by stimulating the left boundary of our simulation domain. We find that, with dx = 0.0225 cm CV increses by 1.3% as we decrease dt from 0.02 to 0.01 ms; if we use dt = 0.02 ms and decrease dx from 0.0225 to 0.015 cm then CV increases by 3.3%; such changes are comparable to those found in earlier studies [15,39].
We introduce conduction inhomogeneities, which we also refer to as obstacles, in the simulation domains of the models described above by making the conductivity constant D = 0 in the region of the obstacle. In most of our studies we use square and squareprism obstacles in two and three dimensions, respectively. When we set D = 0 we decouple the cells inside the obstacle from those outside it. Furthermore, we use Neumann (i.e., no-flux) boundary conditions on the boundaries of the obstacle; we have checked in representative cases that, even if we do not impose Neumann boundary conditions on the obstacle boundaries, our results are not changed qualitatively.
We insert ionic inhomogeneities into our simulation domains by changing the values of the maximal conductances of Ca 2+ channels, in the region of the inhomogeneity, for the LRI, RPB, and TNNP models. Since the Panfilov model does not account explicitly for Ca 2+ ion channels, we mimic ionic inhomogeneities here by altering the value of the parameter e 1 in the region of the ionic inhomogeneity. In most of our studies we use square ionic inhomogeneities in two dimensions.

Results
All the models described above can support spiral waves in a homogeneous simulation domain if we use suitable initial conditions. We begin the description of our results with an overview of such homogeneous simulations for the TNNP model; our work here complements that of Ref. [20]. Spiral waves in homogeneous simulation domains in Panfilov, LRI and RPB models are discussed in Section 3 of the Supplementary Material S1. We then extend our study to simulations with either (a) conduction inhomogeneities or (b) ionic inhomogeneities. This is followed by a discussion of our results on some schemes for the suppression of spiral-wave turbulence in these models; these suppression schemes are the numerical analogs of defibrillation. We consider the efficacy of a few defibrillation schemes in both homogeneous domains and in the presence of the conduction and ionic inhomogeneities described above.

Spiral waves in homogeneous domains
Given the diffusive coupling between cardiac cells, an action potential, generated in one cell, can excite neighboring cells and thus spread out as an expanding wave. However, since the initial excitation is followed by a refractory period, a second wave cannot follow the first one immediately. Each such wave leaves in its wake a nonexcitable region, so the next wave can follow it only at a distance determined by the product of the refractory period and the wave-conduction velocity [40,41]. When two plane waves collide in such a medium, they cannot pass through each other because they have non-excitable wakes. Hence a collision leads to the annihilation of these colliding waves.
Furthermore, the conduction velocities of these waves depend on the curvature of the wavefront [42]: A concave wave moves faster than a rectilinear wave, which in turn travels faster than a convex wave. Any deviations from a planar wave front are, therefore, amplified or attenuated depending on the curvature of the deviation; and they eventually lead to the formation of rotating spiral waves. Above a critical curvature of the wavefront, the current flux from the wavefront is not sufficient to excite the medium around it. This failure of wave conduction then leads to a break up of the wave. These wave fragments move around in the domain, regenerate themselves by using excitable regions, avoid regions that are still in a refractory state, and the parts of these fragments that collide annihilate one another. In a sufficiently large excitable medium this activity of wave fragments can lead to complex spatiotemporal dynamics. The resulting spiral-turbulence state, an instance of spatiotemporal chaos, is characterised by many positive Lyapunov exponents [32].
We show below how such spiral waves can be generated, and how they break up subsequently, in representative, two-dimensional simulations for the TNNP model in homogeneous domains. (Similar simulations for Panfilov, LRI and RPB models are given in Section 3 of the Supplementary Material S1.) Initial conditions have to be chosen carefully to obtain spiral waves; we describe these below. And we use these initial conditions in subsequent Sections that are devoted to our studies of the interactions of such spiral waves with inhomogeneities.
Spiral waves in the TNNP model. We obtain spiral waves in the TNNP model by injecting a plane wave into the domain via a stimulation current of 150 mA/cm 2 for 2 ms at the left boundary. As this plane wave moves towards the right boundary and 270 ms after the first stimulus, we apply a second stimulus of 450 mA/cm 2 along a line behind this wave but parallel to it [14,15] (x = 290, 1#y#250) for 10 ms. As the first wave moves further towards the right, the free end of the new stimulus is able to move into the area behind the first wave; a hook-like proto spiral appears at this free end. We now change the conductivity D from 0.00154 to 0.000385 cm 2 /ms between 304 ms to 524 ms; this yields the fully developed spiral wave. We then reset the conductivity to its original value after 524 ms. At the moment the first plane wave is initiated the currents and gating variables are initialised as follows: the gating variables are given in Supplementary Material S1 and the currents are calculated for these values of the gating variables and the resting value of V which is 286.2 mV. We show the initiation of a spiral wave in the TNNP model in Section 3D of the Supplementary Material S1. The procedure described above results in the spiral wave that is shown via the sequence of pseudocolor plots for the transmembrane potential V (Fig. 1) and the currents I Na , I CaL , I to , I Ks , I Kr , I K1 , I NaCa , I NaK , I pCa , I pK , I bNa , and I bCa (Fig. 2); the states shown in these figures are used as initial conditions for our subsequent simulations of the TNNP model with and without inhomogeneities. In the absence of inhomogeneities such an initial condition leads to a spiral wave as shown in the illustrative pseudocolor plot of V in Fig. 1.
A comparison of spiral waves in different models. From our studies of spiral waves in the four models described above, we see that many qualitative features of spiral-wave dynamics are the same in the Panfilov, LRI, RPB, and TNNP models. However, there are important differences, some qualitative and the others quantitative. For example, the Panfilov model cannot address directly any questions regarding currents in ion channels since it does not follow their evolution but only considers one slow recovery variable g. The LRI, RPB, and TNNP models do give spatiotemporal information for several ion-channel currents as shown, for a representative case, in Fig. 2. At any given time, the qualitative form of the spatial organization of these currents can be surmised from the spatial distribution of the transmembrane potential V and the dependence of these currents on V at the level of a single cell. To illustrate this we show for the TNNP model, in Fig. 3, the temporal evolution of the currents from a single-cell simulation (i.e., without the diffusion term in the TNNP equations). For instance, at the single-cell level, the sodium current I Na is substantial only at the beginning of the action potential (Fig. 3); the most prominent parts of the spiral waves in pseudocolor plots of V appear in regions of the simulation domain where, locally, V assumes a value close to the sharp peak in the single-cell action potential; thus pseudocolor plots of the sodium current I Na (Fig. 2) show significant structure only in narrow strips that follow closely the prominent parts of the spiral waves in pseudocolor plots of V (Fig. 2). By contrast, the calcium current I CaL is significant in the plateau regime of the action potential (Fig. 3); thus pseudocolor plots of I CaL (Fig. 2) show structure in most parts of the simulation domain, but the underlying spiral wave in V is still discernible. The potassium current I K1 is substantial in the repolarization regime of the action potential ( Fig. 3), so we should expect pseudocolor plots of I K1 to have significant structure along the back of this wave, where repolarization occurs; this expectation is borne out as can be seen from Fig. 2. Similar considerations can be used to rationalize, qualitatively, the remaining pseudocolor plots for other currents in the LRI, RPB, and TNNP models; representative plots are given in the Supplementary Material S1.
It is useful to contrast pseudocolor plots of V for the Panfilov, LRI, RPB, and TNNP models. The qualitative features of these plots are the same but they differ in detail [43]. Such differences arises from the differences in the single-cell dynamics of these models, e.g., the action-potential duration, the refractory period, the shape of the repolarization part of the action potential, etc. From the representatives pseudocolor plots of the four models ( Fig. 4), we see that spiral waves in the Panfilov and TNNP models appear more sharp in these plots than their counterparts in the LRI and RPB models. This can be understood qualitatively by comparing the single-cell action potentials for these four models. Figure 5 gives such a comparison, which shows clearly that the repolarization in the Panfilov and TNNP models is sharper and more rapid than in the LRI and RPB models. The sharpness of the spiral waves in the former two models and their relatively diffuse character in the latter two models is related to these differences in repolarization.

Conduction inhomogeneities
We have elucidated the effects of conduction inhomogeneities in the Panfilov and LRI models in Ref. [28]. We begin with a brief recapitulation of these results and then present new ones for the RPB and TNNP models with obstacles. Extensions to the case of two conduction inhomogeneities and their interactions with spiral waves are discussed elsewhere [44,45].
We first examine the dependence of spiral-wave dynamics on the size of an obstacle by fixing its position and changing its size (cf., Ikeda et al.  mm; here we vary l from 2 mm to 80 mm in steps of D = 1 mm. Thus, as l increases, we see a clear transition from ST to RS, with these two states separated by a state Q with no spirals. Henceforth, when we specify the position (x,y) of a square inhomogeneity, we will mean that the bottom-left corner of this inhomogeneity is placed at the point (x,y).
Furthermore, we find that the final state of the system depends on where the obstacle is placed with respect to the tip of the initial spiral wavefront. Even a small obstacle placed near this tip [e.g., an obstacle with l = 10 mm at (100 mm, 100 mm)] can prevent the spiral from breaking up; but a bigger obstacle, placed far away from the tip [e.g., an obstacle with l = 75 mm at (125 mm, 50 mm)], does not affect the spiral. In Ref. [28] we have explored in detail, for the Panfilov and LRI models, how the final state of the system depends sensitively on the position of the obstacle. We have found, in particular, that, if the spiral wave breaks up and yields a spatiotemporally chaotic state in the absence of any obstacles in the medium, then the introduction of an obstacle can lead to one of the following three outcomes: (a) spiral turbulence (ST) can persist; (b) ST can be replaced by a single rotating spiral wave (RS) anchored to the obstacle; (c) ST can give way to a quiescent state (Q) that occurs when all spiral waves move towards and are absorbed by the boundaries.
We show here that this sensitive dependence of spiral-wave dynamics on the position of an obstacle also occurs in the TNNP models. (Similar results for the RPB model are given in Section 4 of the Supplementary Material S1.) We have carried out a systematic study of spiral-wave dynamics in the TNNP model in the presence of an obstacle with the initial conditions specified above, namely, one spiral wave [Figs. 1(A)-(D)] that would continue rotating if the obstacle were not present. In the presence of an obstacle this rotating-spiral (RS) state can be replaced by one of the following possibilities: (a) the spiral wave can continue to rotate, without being anchored to the obstacle, and eventually break down to yield the state ST with spiral turbulence as shown, e.g., in Fig. 6; (b) the tip of the spiral wave can get anchored to the obstacle to give the state RS in which the anchored spiral rotates around the obstacle as shown, e.g., in Fig. 7; (c) all spiral waves can be absorbed by the boundaries so that the system evolves into the quiescent state Q as shown, e.g., in Fig. 8.
Specifically, if we use a square obstacle of side l = 22.5 mm at (54 mm, 22.5 mm), the system is in the state ST as shown in Figs. 6(A)-(C). In this state the time series of the transmembrane potential V (x,y,t), taken from the representative point (90 mm, 90 mm) and shown in Fig. 6(D), clearly displays nonperiodic, chaotic behavior. The time between successive spikes in such a time series, i.e., the interbeat interval (IBI), is plotted versus the integers n, which label successive spikes, in Fig. 6(E); this also shows the chaotic nature of the state ST. We now place the obstacle at (63 mm, 22.5 mm); in this case the spirals get eliminated completely and we get the quiescent state shown via pseudocolor plots of V in Figs. 8(A)-4(C). The time series of V (x,y,t), taken from (x = 90 mm, y = 90 mm) and depicted in Fig. 8(D), clearly shows that we obtain a quiescent state Q with no spirals. Plots of the IBI and the power spectrum are not shown since V just goes to zero after an initial period of transients. The durations for which the transients last, say in Fig. 8(D), vary greatly depending on the position of the obstacle relative to the spiral tip. The sensitive dependence of spiral-wave dynamics on the position of an inhomogeneity is also obtained if we use obstacles that do not have a square shape. We show this   explicitly for a circular obstacle in the Supplementary Material S1 and for obstacles of other shapes in Ref. [45].
This sensitive dependence of the dynamics of spiral waves on the position of an obstacle has been investigated systematically, especially for the Panfilov model, in our earlier work [28]. We have, in particular, obtained a stability diagram for this model as follows: We divide the simulation domain into small squares of side l p = 10 mm. We then carry out a sequence of simulations. In each one of these simulations the square obstacle (slightly larger than the squares into which the simulation domain has been divided) is placed so that its bottom-left corner coincides with the bottom-left corner of one of the small squares; we now study the dynamics of spiral waves, with the initial condition described above, and determine whether the system reaches the ST, RS, or Q state. Our stability diagram, given in Ref. [28], depicts the simulation domain covered with the small squares mentioned above. The color of each small square indicates the final state of the system when the position of the bottom-left corner of the obstacle coincides with that of the small square: red indicates spiral turbulence (ST), blue a rotating spiral (RS) anchored at the obstacle, and green a quiescent state (Q) with no spirals. This stability diagram shows that RS occurs typically when the obstacle lies near the middle of the simulation domain whereas ST occurs when the obstacle lies near the boundary of the simulation domain; regions of Q occur in a few places along the boundary between regions of ST and RS. We have found that this boundary is very complicated; by zooming in on it, we have provided good numerical evidence that suggests that this boundary has a fractal-type character, which leads to the sensitive dependence of the final state of the system on the position of the obstacle. Even if we change the position of the obstacle slightly (say by .0.5 mm), the final state of the system can change from ST to RS or Q. We have suggested [28] that this fractal-type boundary between the ST and RS regions in our stability diagram is a manifestation of an underlying fractal basin boundary between the domains of attraction of the ST, RS, and Q states in the phase space of the infinite-dimensional dynamical system, i.e., the Panfilov-model partial differential equations; such a basin boundary is not easy to determine for an infinitedimensional system but its signatures can be found in the sort of sensitive dependence on parameters, such as the position of the inhomogeneity, that we have elucidated above.
Similar stability diagrams can be obtained, in principle, for LRI, RPB, and TNNP models; however, as the complexity of the models increases so does the difficulty of obtaining a stability diagram. Though we have not found complete stability diagrams for these detailed ionic models, the representative studies that we have carried out indicate that their stability diagrams are qualitatively similar to that of the Panfilov model, which we have given in Ref. [28]. Here we restrict ourselves to parts of such stability diagrams for the TNNP model for the initial conditions described above. (Similar diagrams for the LRI and RPB models are given in Section 4 of the Supplementary Material S1.) For the TNNP model we present a partial stability diagram for G CaL = 0.000044 in Fig. 9, with a square obstacle of side 27 mm, and all other parameters as specified in the figure captions and in the Supplementary Material S1. (Another partial stability diagram for the TNNP model with G pCa = 3.825 is given in Section 4 of the Supplementary Material S1). These partial stability diagrams suggest that the boundaries between ST, RS, and Q states in the TNNP, LRI, and RPB models are as complicated as in the simple Panfilov model. Thus, as we have stated earlier, spiral-wave dynamics in all these models depends very sensitively on the position of a conduction inhomogeneity.
Ionic Inhomogeneities. Apart from obstacles that arise from inhomogeneities in the inter-cellular coupling, cardiac tissue can contain other types of inhomogeneities that originate from changes in single-cell properties, caused say by changes in the chemical environment or metabolic modifications [21,46]. We refer to a collection of such cells, with slightly modified properties like conductances of ion channels, as ionic inhomogeneities; they are different from the obstacles discussed so far in as much as the spiral waves can enter the region spanned by an ionic inhomogeneity. However, such inhomogeneities do affect the dynamics of spiral waves through cardiac tissue [47,48]. In this subsection we investigate spiral-wave dynamics in the presence of ionic inhomogeneities in the four models we have discussed above. For the simple Panfilov model, which does not account for ion channels explicitly, we mimic ionic inhomogeneities via modifications of the inverse time e 1 . We insert such inhomogeneities in LRI, RPB, and TNNP models by considering spatial variations in conductances of calcium ion channels. is one of the recovery time constants [12]. As e 1 increases the absolute refractory period of the action potential decreases. In turn this decreases the pitch of the spiral wave (cf. Fig. 3 in Ref. [35]). Thus by introducing inhomogeneities in e 1 we can investigate spiral-wave dynamics in the presence of time-scale inhomogeneities. If the square simulation domain, of linear size L = 200 mm, is homogeneous, then with e 1 .0.03 we obtain a single, periodically rotating spiral wave but, as it decreases, for instance, if e 1 =0.02, the tip of this rotating spiral starts meandering so that the temporal evolution of the system is quasiperiodic. At even lower values, say at e 1 = 0.01 that we have used above, we see spatiotemporal chaos. These behaviors are shown in the illustrative pseudocolor plots of V in Section 5 of the Supplementary Material S1.
We now introduce a square inhomogeneity in e 1 in the Panfilov model (all other parameters are uniform over the simulation domain): e 1 is assigned the value e in 1 inside a square region; and outside this square it has the value e out 1 . Different choices of e in 1 and e out 1 lead to interesting spiral-wave dynamics. For example, with a square patch of side 40 mm, e in 1~0 :02 and e out 1~0 :01, we obtain spatiotemporal chaos for most positions of this inhomogeneity; but for certain critical positions of this inhomogeneity all spiral waves are completely eliminated; e.g., when the inhomogeneity is at (x = 130 mm, y = 80 mm), spiral waves move towards the boundaries of the simulation domain where they are eventually absorbed. For yet other positions spatiotemporal chaos is obtained outside the inhomogeneity but inside it the spiral wave shows a quasiperiodic temporal evolution. Representative plots are given in Section 5 of the Supplementary Material S1.
If, instead, e in 1~0 :01 and e out 1~0 :02 or 0.03, spiral-wave break up occurs inside the inhomogeneity but it coexists with unbroken  periodically rotating spiral waves outside it (see Section 5 of the Supplementary Material S1) as noted previously by Xie, et al. [31]. However, even in this case, in certain positions such an inhomogeneity anchors a single rotating spiral wave (Fig. 10) as we have seen above, and in Ref. [28], with conduction inhomogeneities; the temporal evolution of V, at a representative point in the simulation domain, is richer than it is with a conduction inhomogeneity: the time series for V can show periodm behavior (we have found cases with 4#m#10) as shown in Figs. 10 A-D for periods m = 5 and m = 4. For example, the plot of the interbeat interval (IBI) versus the beat number n in Fig. 10 D jumps periodically between 167, 280, 244 ms, and 187 ms, i.e., the time series for V displays a period-4 cycle. Such period-m behavior has been reported earlier in experiments [49]; in these experiments it is attributed to the interplay of an anchored spiral wave around very small, but nearby, conduction inhomogeneities.
It is natural to ask whether the rich spatiotemporal dynamics of spiral waves in the Panfilov model with inhomogeneities in e 1 , which lead to inhomogeneities in local times scales since e {1 1 has the dimensions of time, have any analogs in the realistic LRI, RPB, and TNNP models. In cardiac tissue similar changes in time scales can occur because of inhomogeneities in ionic conductances. For example, by changing the conductances G si and G k in the LRI model we can modify the action potential and the time scales associated with it, such as its duration and the extent of the plateau region, which in turn affect spiral-wave dynamics in much the same way as alterations of e 1 do in the Panfilov model. In particular, we can effect transitions from periodic to quasiperiodic rotating spiral waves or from quasiperiodic rotating spiral waves to broken waves with spatiotemporal chaos by changing these conductances; moreover, inhomogeneities in these conductances lead to spatiotemporal patterns in LRI, RPB, and TNNP models that are reminiscent of those we have described above for the Panfilov model with inhomogeneities in e 1 . We show this in detail below by investigating the effects of changes of the maximal calcium and potassium conductances, G si and GK, respectively, in the LRI model, of Gsi in the RPB model, and of the maximal conductance GCaL for the L-type calcium current in the TNNP model. Illustrative simulations for the LRI and RPB models with ionic inhomogeneities are given in Section 5 of the Supplementary Material S1. Here we present our results for ionic inhomogeneities in the TNNP model.
To investigate ionic inhomogeneities in the TNNP model we consider the calcium conductance G CaL that governs the ionic current I CaL (in the Supplementary Material S1 see Section 2 C and Fig. 30 that shows how the action potential (AP) is modified at the single-cell level as we lower G CaL from 0.000175, the maximal channel conductance, to 0.00011, then to 0.00005, and finally to 0 [i.e., I CaL channel block]). We now study the TNNP model in a square simulation domain of side 13.5 cm with the initial condition of Fig. 1 (A). As we decrease G CaL the spiral wave breaks up because the slope of the APD restitution curve steepens and eventually exceeds 1: This break up is shown in the pseudocolor plots of V, at t = 3. and 0.00005, respectively. The discrete lines in the power spectrum of Figs. 11 D can be indexed by one fundamental frequency .8.25 Hz and integer multiples thereof; this is a signature of the periodic rotation of a single rotating spiral wave. The multiple strong peaks and the broad-band background in the power spectra of Figs. 11 E and F are indicative of quasiperiodic (with three fundamental frequencies v 1 .8.25 Hz, v 2 .9 Hz, and v 3 .9.5 Hz) and chaotic states, the latter associated with the break up of spiral waves. We get similar results if we increase the plateau Ca 2+ conductance G pCa instead of changing the L-type Ca 2+ conductance.
We now insert a square G CaL inhomogeneity of side 33.75 mm in a square simulation domain of side 135 mm and G out CaL~0 :000175 (maximal value) and G in CaL~0 :00003, which is approximately one-sixth of its maximal value. When this inhomogeneity is placed at (22.5 mm, 22.5 mm), we observe quasiperiodic behaviors both inside and outside it: Time series for V are recorded from representative points inside and outside of the inhomogeniety, namely, (33.75 mm, 33.75 mm) and (11.25 mm, 11.25 mm), respectively. From these time series we obtain the plots of the IBI and power spectra shown in Fig. 12; we find, in particular, that the main peaks can be indexed as n 1 v 1 +n 2 v 2 +n 3 v 3 with n 1 , n 2 , and n 3 integers and v 1 .4 Hz, v 2 .6.25 Hz, and v 3 .10.5 Hz (inside the inhomogeneity), and v 1 .2.25 Hz, v 2 .4.25 Hz, and v 3 .6.25 Hz (outside the inhomogeneity). Since these frequencies are not related to each other by simple rational numbers we conclude that the spiral wave rotates quasiperiodically both inside and outside the inhomogeneity. In some cases we observe that the inhomogeneity does not have a significant qualitative effect on the dynamics of spiral waves; e.g., when the obstacle is at (45 mm, 45 mm), the position of the spiral tip shifts towards the bottom-left corner of the simulation domain but we still have a state with a single rotating spiral wave whose arms pass through the inhomogeneity. Like conduction inhomogeneities, ionic inhomogeneity can also remove spirals from the medium to leave the system in a quiescent state, e.g., when our G CaL ionic inhomogeneity is at (45 mm, 22.5 mm). (See Fig. 33 in the Supplementary Material S1.)

Elimination of Spiral Turbulence
As we have mentioned above, there is growing consensus that the breakup of spiral waves of electrical activation in ventricular tissue leads to ventricular fibrillation (VF). In the usual clinical treatment of VF electrical stimuli are applied to the affected heart. This is believed to reset all irregular waves in the ventricular tissue leaving it ready to receive the regular sinus rhythm [50]; thus, if These plots indicate that, as G CaL decreases, the system goes from a state with a single rotating spiral wave to the spiral-turbulence state. doi:10.1371/journal.pone.0004738.g011 the electrical stimulus is strong enough, it can arrest VF and restore the sinus rhythm. Initially 60 Hz AC was used clinically to defibrillate transthrorasically [51] but this was later discontinued because of several reasons including the high energy requirement, the possible induction of atrial fibrillation, the prolonged muscle contraction, the risk of an electrical shock to the operator, and the size of the device [50]. Clinically available defibrillation techniques still apply massive electrical shocks to the heart; this can damage the heart muscle. The success rate of such techniques is not quite satisfactory [52]. Furthermore, scar tissues can be created during the process of such defibrillation; these can make the patient vulnerable to further arrhythmias and also act as conduction inhomogeneities that we have investigated via numerical simulations in the Section on ''Results''. Hence there is a great need for developing low-amplitude defibrillation schemes; this must be based on an understanding of the spatiotemporal behavior of activation waves during VF. We begin with a brief overview of some techniques that have been proposed for the elimination of spiral-wave turbulence in models for cardiac tissue. We first examine their efficacy in the homogeneous case; in the next Section we study the relative merits of two low-amplitude defibrillation schemes in the presence of inhomogeneities in our simulation domains. In the context of our numerical simulations of the Panfilov, LRI, RPB, and TNNP models we will use the term defibrillation to mean the elimination of spiral waves or broken spiral waves from the simulation domain.
Early attempts at controlling spiral-wave turbulence in models for cardiac tissue focused on applying well-known techniques of controlling low-dimensional chaos, which are based on the principle that, if a system's trajectory in state space comes very close to an unstable fixed point, it stays in its vicinity for a small duration of time. Different methods have been proposed to drive chaotic trajectories towards such an unstable fixed point, so that the system can stay in the neighborhood of the fixed point for the duration of the control stimulus [53]. But, as we have seen above, the mathematical analog of VF is a state with spiral-wave turbulence. This state displays spatiotemporally chaos, and is, therefore, intrinsically associated with a high-dimensional attractor, so it cannot be controlled by algorithms that have been designed to control chaos in low-dimensional systems.
Review of Previous Work. Biktashev and Holden [54] have proposed a method for controlling spiral-wave turbulence by producing a directed movement of a rigidly rotating spiral wave away from the medium by using resonant stimulation. They find that small-amplitude, spatially uniform, repeated stimuli can be used to produce a directed movement of the spiral wave, if the period of stimulation is equal to the period of its rotation. This directed movement eventually pushes this wave out of the simulation domain [54]. However, this method can only be used before the onset of spiral-wave turbulence.
Osipov and Collins [55] have suggested another scheme that is based on the observation that the dynamics of excitable media can be modelled by fast and slow variables, e.g, V and g in the Panfilov model. They control the slow variable by applying a weak impulse on the whole medium. This eventually changes the velocities of the front and back of the wave. The propagation of the wave front and wave back with different velocities leads to a shrinkage or expansion of the pulse width. If the amplitude and duration of the impulse are sufficiently large, then the propagating pulse collapses and disappears. Unfortunately such control of the slow variable over the whole medium can be achieved only by pharmaceutical means and not by the application of electrical pulses.
Rappel, Fenton, and Karma [56] have proposed another method based on the application of a small control current at a finite number of equally spaced ''controlled cells'' in a tissue, by using a coarse lattice of electrodes with a lattice spacing of about 1 cm. This method has been demonstrated to prevent one spiral from breaking up. Unfortunately this method fails in the fully developed spiral-wave turbulence state with broken spirals [32].
To suppress a spatiotemporally chaotic state with broken spiral waves, Sinha, Pande, and Pandit [32] have proposed a scheme based on the observation that spiral turbulence does not persist in the hearts of small mammals, if it can at all be initiated [57]. We will use this scheme below, so we describe it in some detail. They have shown that spiral-wave turbulence is a long-lived transient [32,35] whose lifetime t L increases rapidly with the linear size L of the simulation domain, e.g., from .850 ms for L = 100 mm to .3200 ms for L = 128 mm in the two-dimensional Panfilov model; for large systems (e.g., L.128 mm in the Panfilov case), t L is sufficiently long so that we obtain a nonequilibrium statistical steady state with spatiotemporal chaos [35]. This might suggest that a global control scheme, such as that of Osipov and Collins [55], is essential. It turns out, however, that a judicious choice of control points on a mesh leads to an efficient scheme for the control of spiral-wave turbulence in such models [32]. We first illustrate the principle of this method for a two-dimensional square domain with side L: This is divided into K 2 smaller blocks by a mesh of line electrodes, and the mesh size is chosen to be small enough that spirals cannot persist for long inside the block of side '~L=K. A voltage or current pulse is applied at all points along the mesh boundaries for a time. This makes the mesh region refractory and so effectively simulates Neumann boundary conditions for any block bounded by the mesh. Thus spiral waves formed inside the block are absorbed at the mesh bounding the block. For example, in the Panfilov model in dimension d = 2, L = 128, and K = 2, a time t c~4 1:2 ms suffices to suppress spiral turbulence; when L = 512, and K = 8, a time t c~7 04 ms is required; electrical pulses of amplitude .60 mA/cm 2 are used on the control mesh; this is much less than in conventional electrical defibrillation which uses pulses of amplitude 1 A/cm 2 . This control algorithm has been extended to suppress spiral turbulence in the two-dimensional Beeler-Reuter and LRI models [35]. A nave extension of this control algorithm to three-dimensions requires a cubic array of sheets, which cannot be implanted easily in ventricular tissue. However, in Ref. [32] it was shown that, even if the control mesh is present only on one of the L6L faces of an L6L6L z domain, the above scheme works with a slight modification: Instead of applying a pulse for a time t c , we have to apply a sequence of n pulses, separated by t ip ; the duration of each pulse is t v . A steady pulse does not control scroll-wave turbulence in three dimensions for L z .2 mm: its propagation in the z direction is impeded once the interior of the simulation domain becomes refractory. However, if we use a sequence of short pulses and if t ip is long enough for the medium to recover its excitability, the control-pulse waves can propagate in the z direction and lead to successful control. We refer the reader to Refs. [32,35] for further details; we give representative threedimensional simulations after our discussion of two-dimensional studies.
Recently Zhang, et al. [33], have proposed another attractive scheme for the control of spiral turbulence in excitable media. In their method spiral waves are driven away by periodic forcing of V at a small number of n6n points in the center of the simulation domain. This generates target waves that eventually drive out the spiral waves if the amplitude C and the frequency v f of the forcing are chosen carefully: For example, for the Panfilov model with d = 2 it is shown in Ref. [33] that spiral turbulence in a square 5006500 simulation domain can be suppressed within 410,000 iterations when one chooses n = 6, v f = 0.82, and C~6. This control scheme is attractive because it employs local forcing, compared to the control scheme of Ref. [32] that uses a spatially extended control mesh. However, as we show below, the local control scheme of Ref. [33] inadvertently generates spiral-wave break up if there are obstacles in the medium.
Summary of Our Results. The control scheme of Ref. [32] is also successful in eliminating spiral turbulence in the realistic LRI, RPB, and TNNP models that account for ion channels. Illustrative simulations for the LRI and RPB models are given in Section 6 of the Supplementary Material S1. Figure 13 presents results from our simulations of the two-dimensional TNNP model; here a control current of 27.75 mA/cm 2 , applied for 20 ms on a mesh that divides our square simulation domain of side 13.5 cm into 16 square cells of side 3.375 cm each, suffices to control spiral turbulence. The pseudocolor plots of V in Fig. 13 give (A) the initial spiral-turbulent state and its subsequent evolution after the initiation of the control pulse at (B) 24 ms, (C) 80 ms, and (D) 280 ms, after which all spiral turbulence disappears.
As we mentioned above, the control scheme of Ref. [32] can be extended to three-dimensions. For a practical implementation of this scheme the control pulse must be applied only on one face of the simulation domain. This suffices when the thickness in the third dimension is small. For example, in the three-dimensional Panfilov model, even if we apply a control pulse on a mesh on one of the L6L faces of the L6L6L z simulation domain, we can suppress broken scroll waves (the three-dimensional analogs of broken spiral waves) in the entire medium, provided that L z ,2 mm. Illustrative simulations are given in Section 6 of the Supplementary Material S1.
For L z .2 mm this control scheme for a three-dimensional domain must be modified as follows: We use a control mesh on one face, but, instead of using a single long pulse, we use a series of short pulses, each of duration t v and separated by an interval t ip ; for suitable choices of t v and t ip , we can effectively suppress spiral turbulence. (A similar scheme, with short pulses separated from each other, also works in two dimensions.) We illustrate this control scheme for the three-dimensional Panfilov model in a 128612864 mm 3

The Control of Spiral-Wave Turbulence in the Presence of Inhomogeneities
In the last Section we have given a short overview of some control schemes that have been used to suppress spiral-wave break up in two-and three-dimensional simulations in some mathematical models for cardiac tissue. Cardiac tissue can have inhomogeneities, such as scar tissue. It is important, therefore, to study whether these control schemes are effective in controlling spiral-wave turbulence in the presence of such inhomogeneities. To the best of our knowledge this has not been investigated systematically so far. We present such an investigation here (for simplicity we restrict ourselves to conduction inhomogeneities). In particular, it is important to ensure that a control scheme does not lead inadvertently to spiral break up in the presence of inhomogeneities.
We begin by studying the control scheme of Zhang et al.
[33] that we have outlined above. This scheme drives away broken spiral waves from the simulation domain by using the target waves that are created by the local periodic forcing. What happens to such target waves when they encounter an obstacle? We examine this for the two-dimensional Panfilov model. In Fig. 15 we present our results with such periodic forcing [Ccos v f t ð Þ with C~6 and v f = 0.82]; we use parameters such that, in the absence of this forcing, a single spiral wave would have been attached to the conduction inhomogeneity ( Fig. 15 A that is similar, e.g., to Fig. 1B in Ref. [28]). The forcing we use generates target waves in the center of the simulation domain, of side L = 25 cm and with a square obstacle of side l = 4 cm placed at (8.5 cm, 8.5 cm). It turns out that these target waves drive away the spiral wave that was anchored to the obstacle in the absence of the forcing. The time evolution of the system with the forcing (Figs. 15 B-D) shows the break up of these target waves as they collide with the obstacle and thus contribute to spiral turbulence in the medium [58]. Had there been no obstacle, this control scheme would have driven away all the broken spiral waves from the domain. However, this does not happen in the presence of an obstacle as shown in Fig. 15; hence this control scheme is unsuitable for controlling spiral-wave turbulence if inhomogeneities are present.
By contrast, the control scheme proposed in Ref. [32] works even in the presence of an inhomogeneity. Illustrative simulations for the Panfilov, LRI, and RPB models are given in Section 6 of the Supplementary Material S1. Here we show how the control scheme of Ref. [32] is also successful in eliminating spiral turbulence in the TNNP model even in the presence of conduction inhomogeneities. In Fig. 16 a control current of 27.75 mA/cm 2 is applied for 20 ms on a mesh that divides our square simulation domain of side 13.5 cm into 16 square cells of side 3.375 cm each. This suffices to control spiral turbulence in the TNNP model, even though there is a square obstacle of side 2.25 cm placed in the simulation domain. (If this obstacle is at (45 mm, 31.5 mm), in the absence of the control pulse, a single rotating spiral wave would have been anchored to this obstacle; if instead the obstacle is at (54 mm, 22.5 mm), as in Fig. 6, spirals would have broken up in the absence of the control.) The pseudocolor plots of V in Fig. 16 give (A) the initial state and its subsequent evolution after the initiation of the control pulse at (B) 24 ms, (C) 80 ms, and (D) 280 ms, after which all spiral turbulence disappears. A similar simulation for the elimination of spiral turbulence in the TNNP model is given in Section 6 of the Supplementary Material S1.
We have shown above how the control scheme of Ref. [32] can be extended to three-dimensions in an L6L6L z simulation domain, if L z ,2 mm. This scheme works even in the presence of an obstacle as is illustrated for representative cases in Fig. 17 with an obstacle of size 2562562 mm 3 placed at (80 mm, 60 mm). In the absence of control pulses, scroll-wave turbulence is obtained. We now apply a control pulse of strength 0.48 on a mesh with square cells, each of side 16 mm, on the bottom face of the simulation domain for 748 ms; the evolution of the states of the system after the initiation of the control are depicted at (B) 220 ms, (C) 440 ms, and (D) 880 ms in Fig. 17. By 880 ms the scroll waves have left the simulation domain and it is completely quiescent.
When L z .2 mm, the control scheme described in the previous paragraph fails just as its counterpart did in the absence of an inhomogeneity. A representative simulation for this case is given in Section 6 of the Supplementary Material S1.

Discussion
We have presented the most extensive numerical study carried out so far of the effects of inhomogeneities on spiral-wave dynamics in mathematical models for cardiac tissue. In particular, we have investigated such dynamics in the Panfilov, LRI, RPB, and TNNP models for homogeneous simulation domains and also in the presence of conduction and ionic inhomogeneities. Furthermore, we have considered two low-amplitude control schemes in detail; these have been designed to eliminate spiralwave turbulence in these models but have not been tested systematically in the presence of inhomogeneities; we carry out such tests here.
One of the principal results of our studies is the confirmation that spiral-wave and scroll-wave dynamics in mathematical models of cardiac tissue depend very sensitively on the positions of conduction or ionic inhomogeneities in the simulation domain. Our results here extend significantly those presented in Ref. [28] for the Panfilov and LRI models. In particular, we have shown that this sensitive dependence on inhomogeneities also holds in realistic ionic models, which account for ion pumps and ion exchangers and also the details of the dynamics of calcium ions; furthermore, these results also hold in three-dimensional simulation domains, as illustrated by our calculations for the threedimensional Panfilov model; and the nature of the inhomogeneity also affects the spatiotemporal dynamics of spiral waves as can be seen by comparing our simulations of conduction inhomogeneities with those for ionic inhomogeneities. As we have seen, in the latter case the transmembrane potential V displays rich and different temporal behaviors inside and outside the ionic inhomogeneity. We believe this sensitive dependence of spiral waves on inhomogeneities in the medium is a reflection of a fractal basin boundary between the domains of attraction of spiral-turbulence (ST), rotating-spiral (RS), and quiescent (Q) states. In a lowdimensional dynamical system it is possible to obtain such a basin boundary by changing initial conditions; in a high-dimensional dynamical system (the partial-differential-equation models for cardiac tissue are infinite dimensional) it is not practical to find such a boundary numerically. We have shown instead, that, by changing parameters in these cardiac-tissue models, such as the positions or natures of inhomogeneities, we can affect the spatiotemporal evolution of spiral waves drastically.
Our studies have practical implications for experimental investigations of spiral-wave dynamics in cardiac tissue. In particular, the studies of Refs. [22][23][24]49] have provided a rich variety of results including complicated temporal patterns in inter-beat intervals [49] for V and the partial or complete elimination of spiral-wave turbulence by conduction inhomogeneities [2]. We have described these briefly in the introduction. Here we would like to note that our in silico simulations of spiralwave dynamics in the Panfilov, LRI, RPB, and TNNP models have allowed us to carry out a much more systematic study of inhomogeneities in these models than is possible in vitro and in vivo. We hope our work will stimulate experiments in this field. It is worth noting that our study yields all the types of rich spatiotemporal behaviors (e.g., for V) that have been observed in a variety of experiments on spiral-wave dynamics in cardiac tissue or cell cultures, if we keep in mind that the states ST, RS, and Q in our simulations are the analogs of VF, VT, and quiescence in such experiments.
Our results, especially those on the elimination of spiral-wave turbulence in the presence of inhomogeneities, should also have important implications for the development of low-amplitude electrical defibrillation schemes, which is a major challenge that lies at the interfaces between nonlinear science, biophysics, and biomedical engineering. One of the lessons of our numerical studies, namely, the sensitive dependence of spiral-wave dynamics on inhomogeneities, implies that low-amplitude defibrillation schemes might well have to be tuned suitably to account for inhomogeneities in cardiac tissue. Furthermore, it would be very interesting to develop the mesh-based control scheme that we have described in the previous section and to see how it might be realised experimentally.
While this paper was being prepared for publication two new suggestions appeared for the elimination of spiral turbulence in models such as the LRI model [59,60]. The first of these [59] uses a square lattice of control points through which a control voltage is swept. Since this uses a spatially extended set of control points, it is successful in the elimination of spiral turbulence at least in a twodimensional domain. The study of Ref. [60] examines a scheme that requires the use of a bidomain model since relatively large control voltages are used. We have concentrated instead on lowamplitude defibrillation schemes that should not require bidomain models; indeed such schemes have been used earlier by several groups [32,33,56,61,62] without bidomain models.
As we have emphasized throughout this paper, one of the principal goals of our study is a qualitative one, namely, the elucidation of the sensitive dependence of spiral-wave dynamics on inhomogeneities in mathematical models of cardiac tissue. We have, therefore, carried out extensive simulations of such dynamics in the Panfilov, LRI, RPB, and TNNP models; but we have not, so far, extended our study to bidomain models [63] and models in which mechanics [64] is also included. We expect our principal qualitative results about inhomogeneities will go through even when such models are considered; this will have to be checked explicitly by subsequent studies.
In our earlier 3D simulations [28] we have studied the effect of inhomogeneities on scroll waves; in this case too we find that scroll-wave dynamics depends sensitively on the position of an inhomogeneity. Here too we discuss some 3D simulations. However, we have concentrated on 2D simulations for two reasons: (a) it is important to have a good understanding of spiral-wave dynamics in the presence of inhomogeneities before embarking on similar, detailed, 3D simulations; (b) in 3D studies it is important to include tissue anisotropy (but again we believe this can be done systematically only after we have a good understanding of the 2D simulations we have presented here). In the same spirit of studying the simplest models first, we have not considered models which include mechanics too; indeed this purely electrical approach has been adopted by several other studies (see, e.g., Refs. [11,12,14,25,[30][31][32][33]) with the understanding that the mechanical system basically follows the electrical activation at the level of a first approximation (see, e.g., Ref. [60]).

Supporting Information
Supplementary Material S1 Spiral-Wave Turbulence and its Control in the Presence of Inhomogeneities in Four Mathematical Models of Cardiac Tissue Found at: doi:10.1371/journal.pone.0004738.s001 (3.99 MB PDF)

Author Contributions
Conceived and designed the experiments: STK RP. Performed the experiments: STK ARN RP. Analyzed the data: STK ARN RP. Contributed reagents/materials/analysis tools: STK ARN RP. Wrote the paper: STK ARN RP.