Spiral-wave dynamics in ionically realistic mathematical models for human ventricular tissue: the effects of periodic deformation

We carry out an extensive numerical study of the dynamics of spiral waves of electrical activation, in the presence of periodic deformation (PD) in two-dimensional simulation domains, in the biophysically realistic mathematical models of human ventricular tissue due to (a) ten-Tusscher and Panfilov (the TP06 model) and (b) ten-Tusscher, Noble, Noble, and Panfilov (the TNNP04 model). We first consider simulations in cable-type domains, in which we calculate the conduction velocity θ and the wavelength λ of a plane wave; we show that PD leads to a periodic, spatial modulation of θ and a temporally periodic modulation of λ; both these modulations depend on the amplitude and frequency of the PD. We then examine three types of initial conditions for both TP06 and TNNP04 models and show that the imposition of PD leads to a rich variety of spatiotemporal patterns in the transmembrane potential including states with a single rotating spiral (RS) wave, a spiral-turbulence (ST) state with a single meandering spiral, an ST state with multiple broken spirals, and a state SA in which all spirals are absorbed at the boundaries of our simulation domain. We find, for both TP06 and TNNP04 models, that spiral-wave dynamics depends sensitively on the amplitude and frequency of PD and the initial condition. We examine how these different types of spiral-wave states can be eliminated in the presence of PD by the application of low-amplitude pulses by square- and rectangular-mesh suppression techniques. We suggest specific experiments that can test the results of our simulations.


INTRODUCTION
Sudden cardiac arrest is caused, in many cases, by cardiac arrhythmias, such as ventricular tachyacardia (VT) and ventricular fibrillation (VF) (Roger et al., 2011(Roger et al., , 2012. Estimates suggest that VF is the main reason for death in 30% of the cases in which heart failure occurs (Zipes and Wellens, 1998;Fogoros, 2011). Thus, the importance of studying such arrhythmias cannot be overemphasized. Such studies must use interdisciplinary approaches because they require inputs from biology, bio-medical engineering, and cardiology, on the one hand, and physics, non-linear dynamics, and numerical methods, on the other; methods from these areas must be used to study the complicated, non-linear, partialdifferential-equation models that have been developed for cardiac tissue. Such equations can show, inter alia, spiral-wave turbulence and spatiotemporal chaos, which is believed to be one of the mathematical analogs of VF. The study we present here combines theoretical ideas from spatiotemporal chaos in extended dynamical systems with extensive direct numerical simulations, to elucidate the effects of periodic deformation (PD) on spiral-wave dynamics in detailed mathematical models for cardiac tissue and to investigate the elimination of such spiral waves, in the presence of PD, by the application of low-amplitude current pulses.
The mechanisms underlying VT and VF are not understood with complete certainty; however, various clinical studies Fogoros, 2011) have suggested that such arrhythmias comprise the abnormal propagation of a wave of electrical activation across the ventricles. Such irregular waves appear in various ways. For example, they can arise because of an infarction scar (De Bakker et al., 1993), which can create an anatomical path that anchors electrical waves; they are also seen around obstacles in vitro (Valderrábano et al., 2000) and in cell-culture experiments (Lim et al., 2006) (anatomically reentry). However, reentry can occur in the absence of an anatomical pathway. For instance, in vitro experiments (Davidenko et al., 1990;Ikeda et al., 1996) have shown that fibrillation can be maintained by functional reentry. However, ectopic activations (Chen et al., 1999;Zimmermann and Kalusche, 2001), because of pulmonary veins or abnormal cells, can also initiate fibrillation. In particular, both experimental (Davidenko et al., 1992;Pertsov et al., 1993;Gray et al., 1998;Jalife et al., 1998) and computational Fenton et al., 2002;Cherry and Fenton, 2008) studies have suggested that VT and VF are, respectively, manifestations of (a) a rotating spiral (RS) or scroll wave or (b) broken spiral or scroll waves that lead to spiral-or scroll-wave turbulence (ST).
Several studies have investigated the transition from RS to ST, both in experiments on cardiac tissue and in computational studies of mathematical models for cardiac tissue; they find that this transition can occur because of (a) a steep, increasing initial segment in the restitution curve, a plot of the action potential duration (APD) versus the diastolic interval (DI) (Koller et al., 1998;Garfinkel et al., 2000;Fenton et al., 2002), (b) a similar steep part in an analogous plot of the conduction velocity θ versus DI (Qu et al., 2000b;Fenton et al., 2002), (c) alternans (Karma, 1994;Koller et al., 1998;Qu et al., 1999;Cherry and Fenton, 2008), and (d) heterogeneities, such as, conduction and ionic inhomogeneities (Xie et al., 1998;Shajahan et al., 2007Shajahan et al., , 2009Majumder et al., 2011a,b). Recently, some groups (Zhang et al., 2004(Zhang et al., , 2006Panfilov et al., 2007;Chen et al., 2008;Weise et al., 2011) have begun to study the effects of the deformation of cardiac tissue on the RS-ST transition; such a transition arises either because of periodic deformation or the stretchactivated current associated with such deformation. These studies have used simple, two-variable mathematical models for electrical activation in such tissue. One of the goals of our study is to investigate spiral-wave dynamics in general, and RS-ST transitions in particular, in a simple mathematical model for periodic deformation (PD) of cardiac tissue (Zhang et al., 2004(Zhang et al., , 2006Chen et al., 2008) that we couple with ionically realistic human-ventricular-tissue mathematical models, namely, (a) the TP06 model, due to ten Tusscher and Panfilov (Ten Tusscher and Panfilov, 2006), or (b) the TNNP04 model, of ten Tusscher, Noble, Noble, and Panfilov (Ten Tusscher et al., 2004). A deformation of cardiac tissue leads to modifications of ion-channel parameters, because of stretch-activated currents, and also a modification of intracellular couplings. Our model for deformation, based on Zhang et al. (2004Zhang et al. ( , 2006 and Chen et al. (2008), is a very simplified one in which the effects of deformation are accounted for only by a temporal modulation of diffusion constants (see section 2), which are related to intracellular couplings; we do not include stretch-activated currents as considered in Panfilov et al. (2005Panfilov et al. ( , 2007 and Weise et al. (2011). However, in spite of this simplified representation of deformation, our study yields important results that have been observed in twovariable models for cardiac tissue both with periodic deformation (PD) (Zhang et al., 2004(Zhang et al., , 2006Chen et al., 2008) or mechanical deformation (Panfilov et al., 2005(Panfilov et al., , 2007Weise et al., 2011); the latter studies include stretch-activated currents. Such stretchactivated currents with mechano-electrical feedback have also been shown to affect electrical activation in anatomically realistic human models (Keldermann et al., 2010;Kuijpers et al., 2011); mechano-electrical feedback can enhance electrical activation as discussed in Thompson et al. (2011) in the context of myofibroblast-myocyte interactions (of course, pharmacological and electrochemical interventions can also enhance electrical activation). On the positive side, our study uses ionically realistic models that have not been employed in such deformation studies so far (but see the recent Weise and Panfilov, 2013). We discuss the principal results of the studies carried out in Zhang et al. (2004Zhang et al. ( , 2006, Panfilov et al. (2005Panfilov et al. ( , 2007, Chen et al. (2008), and Weise et al. (2011), in section 4, where we compare their findings with ours.
We also investigate the efficacy of a low-amplitude suppression scheme, developed for the suppression of spiral-wave turbulence in 2D models for cardiac tissue (Sinha et al., 2001;Pandit et al., 2002;Shajahan et al., 2009;Majumder et al., 2011a) in the absence of PD. Our study yields several interesting results that we summarize below before we discuss them in detail. We find first that PD leads to a periodic, spatial modulation of the conduction velocity θ and a temporally periodic modulation of the wavelength λ of a plane wave. We then use three different parameter sets, for both TP06 and TNNP04 models, to obtain three different prototypical spiral configurations, which we use as the initial conditions IC1, IC2, and IC3 (see section 2). We find, for the TP06 model, that spiral-wave dynamics depends sensitively on PD and the initial condition. (For similar studies of the sensitive dependence of spiral-wave dynamics on inhomogeneities, see Shajahan et al., 2007Shajahan et al., , 2009.) The initial condition IC1 can lead to (a) an RS state with n-cycle temporal evolution (here n is a positive integer), (b) rotating-spiral states with quasi-periodic (QP) temporal evolution, (c) a state with a single meandering spiral MS, which displays spatiotemporal chaos, (d) an ST state, with multiple broken spirals, and (e) a quiescent state SA, in which all spirals are absorbed; the initial condition IC2, with PD, can lead either to (a) an ST state, with multiple spirals, or (b) an SA state, with no spirals; and for IC3, it can be driven into (a) an ST state, with a single meandering spiral, (b) an ST state, with multiple spirals, and (c) the state SA. For all these initial conditions, precisely which one of these states is obtained depends on the amplitudes A x and A y and the frequencies f x and f y of the PD in the x and y directions. Spiral-wave dynamics in the TNNP04 model, with PD, also shows sensitive dependence on PD and the initial condition. This sensitive dependence on parameters is a hallmark of extended dynamical systems that show spatiotemporal chaos (Shajahan et al., 2007(Shajahan et al., , 2009. We also study, in the presence of PD, the efficacy of a low-amplitude suppression scheme (Sinha et al., 2001;Shajahan et al., 2009) that has been suggested, hitherto only without PD, for the suppression of spiral-wave turbulence, via low-amplitude current pulses applied on a square mesh, in mathematical models for cardiac tissue. Furthermore, we develop line-mesh and rectangular-mesh variants of this suppression scheme. The latter suppresses spiral turbulence in all cases we consider.

METHODS
The electrical activation of the transmembrane potential V m of cardiac tissue is often modeled by a reaction-diffusion-type equation, where C m is the membrane capacitance density, I ion is the sum of all the ionic currents that cross the cell membrane, and D x and D y are, respectively, the diffusion coefficients along x and y directions; such diffusion terms are related to gap junctions (Ten Tusscher et al., 2004;Ten Tusscher and Panfilov, 2006), which are networks of protein channels that allow the passage of ions from cell to cell. We use two biophysically realistic ionic models for human cardiac myocytes: (a) the ten Tusscher and Panfilov model (the TP06 model) (Ten Tusscher and Panfilov, 2006) and (b) the ten Tusscher, Noble, Noble, and Panfilov model (the TNNP04 model) (Ten Tusscher et al., 2004). It is important to check the compatibility of the models that we consider in our study with the standard version of the original models (Ten Tusscher et al., 2004;Ten Tusscher and Panfilov, 2006). One way to check for such compatibility is by comparing the AP and its morphological properties ( Figure S1 and Table S1, Supplementary Material S1) with those in our studies; we have checked this explicitly, as we discuss in detail in the Supplementary Material S1. We follow the method suggested in Zhang et al. (2004Zhang et al. ( , 2006 and Chen et al. (2008) for the introduction of PD into a mathematical model for cardiac tissue. In particular, we note that any point x = (x, y) in the medium changes to x (t) = (x (t), y (t)) with and A y (t) = A y cos (2π f y t). By substituting Equation (2) into Equation (1), we obtain a comparison of Equations (1, 3) shows that Equation (3) can be rewritten as with D x (t) = D x (1 + A x (t)) −2 and D y (t) = D y 1 + A y (t) −2 .
In our numerical simulations, we use 2D square domains with 1024 × 1024 grid points and lattice spacings δx = δy = 0.25 mm for both TP06 and TNNP04 models, so the sides of our square simulation domains are L = 256 mm in the absence of PD. We use a forward-Euler method for time evolution, with a time step δt = 0.02 ms, a five-point stencil for the Laplacian, and no-flux (Neumann) boundary conditions. We set the diffusion coefficients D x = D y = D = 0.00154 cm 2 /ms (Ten Tusscher et al., 2004;Ten Tusscher and Panfilov, 2006) for both the TP06 and the TNNP04 models for our numerical investigations. Other parameters for our calculations are given in Tables 1, 2 and an examination of the numerical stability of our numerical scheme is given in the Supplementary Material S1.
To examine the spatiotemporal evolution of the electrical signal of the transmembrane potential, we obtain the local time series of V m (x, y, t), from a representative point (x = 125 mm, y = 125 mm) (shown by an asterisk in all pseudocolor plots of V m ). To obtain the plots of the inter-beat interval (IBI), we use this local time series with 4 × 10 5 data points; the IBI is the time interval between two successive beats in the time series of the signal V m . For the power spectra E(ω) we use the local time series with 2 × 10 5 data points after the initial 10 5 data points have been removed to eliminate transients. We present the spatiotemporal evolution of V m by a series of (Videos S1-S9) of its pseudocolor plots; all these videos use 10 frames per second and each frame is separated from the succeeding frame by 8 ms.
We often have to track the trajectory of the tip of a spiral wave in a 2D simulation domain. The tip of such a spiral wave is normally defined as the point where the excitation wave front Here, σ f is the scale factor of the time constant τ f (Ten Tusscher et al., 2004;Ten Tusscher and Panfilov, 2006). and repolarization wave back meet; this point can be found by a variety of methods (Barkley et al., 1990;Fenton and Karma, 1998;Fenton et al., 2002;Otani, 2002;Gray et al., 2009;Nayak et al., 2013). We use the tip-tracking algorithm of Nayak et al. (2013) that locates the tip position by monitoring I Na , the sodium current. Pseudocolor plots of I Na show a fine line along the arm of a spiral wave ( Figure 2A in Shajahan et al., 2009); this line terminates in the spiral tip and can, therefore, be used to obtain the spatiotemporal evolution of this tip.
In Figures

RESULTS
We begin by exploring the effects of PD both on plane-wave propagation and on spiral-wave dynamics; here we vary the oscillation amplitude and the frequency in the ranges 0 A x , A y 0.5 and 0 Hz f 7.0 Hz; the deformation amplitudes we use are comparable to those in other computational (Zhang et al., 2004(Zhang et al., , 2006Weise et al., 2011) andexperimental (McCULLOCH et al., 1987;Noble, 2002) studies; to set the scale of frequencies, we note that the frequency of rotation of a single spiral wave is 4.75 Hz for the TP06 model and 3.75 Hz for the TNNP04 model (see section 3.2). We then study the effects of PD on the suppression scheme of Pandit et al. (2002) and Shajahan et al. (2009).

PLANE-WAVE DYNAMICS IN A CABLE
We study plane-wave propagation in a thin, cable-type simulation domain, with 16 × 4096 grid points, i.e., L x = 4 mm and L y = 1024 mm. We inject a stimulus of strength I stimulus = 150 pA/pF at the left end of the cable for 3 ms and then study the effects of PD on the plane wave that propagates through this cable; in particular, we measure the conduction velocity θ and wavelength λ FIGURE 1 | Pseudocolor plots of the transmembrane potential V m for the TP06 model illustrating plane-wave propagation in a cable-type domain, with PD along the axial-direction of the cable, and the parameter sets given in Table 1. The Video S1 comprises 21 animations that show the spatiotemporal evolution of these plane waves.
of the propagating wave in the cable. We find that θ 70.6 cm/s and λ 21.6 cm for the TP06 model, and θ 67.8 cm/s and λ 18.9 cm for the TNNP04 model, in the absence of PD. As suggested in Clayton and Panfilov (2008), Shajahan et al. (2009), andTen Tusscher et al. (2004), it is useful to test the accuracy of the numerical scheme by varying both the time and space steps that we use for integration. We illustrate this for the TP06 model by measuring θ for a plane wave, which is injected into the medium by stimulating the left boundary of our simulation domain. We find that, with δx = 0.025 cm, θ increases by 1.6% as we decrease δt from 0.02 to 0.01 ms; if we use δt = 0.02 ms and decrease δx from 0.025 to 0.015 cm then θ increases by 4.7%; such changes are comparable to those found in earlier studies (Ten Tusscher et al., 2004;Shajahan et al., 2009).
In Figure 1(a00-a20) we show, at time t = 600 ms, when PD is applied along the axial direction of the cable, pseudocolor plots of the transmembrane potential V m for the TP06 model with PD along the axial direction of the cable, and the parameter sets given in Table 1. The Video S1 comprises 21 animations that show the spatiotemporal evolution of the plane waves in Figure 1(a00-a20); these animations and Figure 1(a00-a20) show that the conduction velocity θ is modulated in space and the wavelength λ is modulated in time because of the PD. Figure 2 illustrates these modulations via plots of θ F and θ B versus x for the conduction velocities of the wave front (Figure 2A) and the wave back ( Figure 2B), respectively; here the subscripts F and B stand for wave front and wave back, respectively; and Figure 2C shows the corresponding plot for λ versus time t; in these plots we use the representative PD parameter values A x = 0.3 and f x = 5.0 Hz for the TP06 model. We calculate the conduction velocities θ F (x) and θ B (x), in the cable-type domain with PD, by recording the positions of the wave front and the wave back at times t and t + t, with t = 2 ms; the wave-front and wave-back conduction velocities, at the point x at time t, where F x and B x are, respectively, the distances traveled by the wave front and wave back in the time interval t. We locate the position of the wave front by finding the value of x at which V m 0 mV; we define the position of the wave back as the point, behind the wave front, at which a secondary action potential can just be initiated by an additional stimulus (this turns out to occur at a value of V m that is 75% of the repolarization phase of the action potential). We obtain the wavelength λ(t) by measuring the distance between the wave front and the wave back at time t.
In Figures 2A,B, the open circles show the values of θ F (x) and θ B (x), respectively, that we obtain by the method described above; the red lines show smooth sinusoidal envelopes, with amplitude 31.2 cm/s and spatial period 14.5 cm, that give the average modulations of these conduction velocities with x. Note that, in the absence of PD, θ 70.6 cm/s (this is shown via a gray, dashed line in Figures 2A,B); therefore, the electrical wave can travel 70.6/f cm in 1/f s; hence, for a given PD frequency f , the spatial period of oscillation of θ F (x) and θ B (x) is 70.6/f cm; the representative plots of Figures 2A,B, in which f = 5 Hz and the period is 70.6/5 = 14.12 cm, are consistent with this estimate. Figure 2C shows that λ is a periodic function of t with a period τ ; we expect that τ = 1/f , where f is the PD frequency; the illustrative plot in Figure 2C, with f = 5 Hz, is consistent with this expectation because τ 202 ms; the gray, dashed line shows the value of λ that we obtain in the absence of PD.
It is useful to study how θ and λ of a plane wave behave, in the presence of PD, when we change the values of the time step and the diffusion coefficients. We find that, in the presence of PD, θ and λ continue to oscillate, as in Figures  The TNNP04-model analogs of the TP06-model Figure 1(a00-a20) are given in Figure S4(a00-a11) in the Supplementary Material S1. Our results for the TNNP04 model are similar to those for the TP06 model.

SPIRAL-WAVE DYNAMICS IN A HOMOGENEOUS DOMAIN
We move now to systematic studies of spiral-wave dynamics in a 2D, square simulation domain with side L = 256 mm, in the presence of PD, for both TP06 and TNNP04 models.
We describe in the Supplementary Material S1 the precise S1, S2 cross-field protocol that we use to obtain spiral waves in our simulation domain. We use three types of spiral-wave initial configurations for our subsequent studies; we refer to these as IC1, IC2, and IC3 initial conditions (see Table 2 for parameter values). In Figures  In Figures 3A-C, we show pseudocolor plots of V m at times t = 0 s, t = 2 s, and t = 4 s, respectively, for the initial condition IC1 in the TP06 model, in the absence of PD; this initial configuration evolves to a state with a rotating spiral (RS) in the medium; the animation (a) in Video S2 shows the spatiotemporal evolution of V m for this case. The local time series of V m (x, y, t), from the representative point (x = 125 mm, y = 125 mm) (the asterisk in Figure 3C), is shown in Figure 3D for 2 s ≤ t ≤ 6 s; a plot of the IBI is given in Figure 3E, which shows that, after initial transients (roughly the first 10 beats), the spiral wave rotates periodically with an average rotation period T 210 ms. In Figure 3F, we plot the power spectrum E(ω), which we have  (x, y, t), from the representative point (x = 125 mm, y = 125 mm) (the asterisk in (C)) for 2 s ≤ t ≤ 6 s; (E) a plot of the IBI, which we obtain from this time series, of length 4 × 10 5 iterations; (F) the power spectrum E(ω), obtained from the local time series of (D), with discrete peaks at the fundamental frequency ω f 4.75 Hz and its harmonics. The spiral-tip trajectory traces a roughly circular path, with radius l c 20 mm, which is shown, for 3.6 s ≤ t ≤ 4 s, by the white line that has been superimposed on the pseudocolor plot of V m in (C); a magnified view of this path is shown in (G). obtained from the local time series of V m mentioned above; discrete peaks in E(ω) appear at the fundamental frequency ω f 4.75 Hz and its harmonics. The periodic nature of the local time series of V m , the flattening of the IBI, and the discrete peaks in E(ω) show that the temporal evolution of the spiral wave is periodic; therefore, the spiral-tip trajectory traces a roughly circular path with radius l c 20 mm; this circular path is shown, for 3.6 s ≤ t ≤ 4 s, by the white line that has been superimposed on the pseudocolor plot of V m in Figure 3C; an expanded version of this path is shown in Figure 3G. The Figures S7, S8 in the Supplementary Material S1 show the analogs of Figure 3 for initial conditions IC2 and IC3, respectively; and the animations (b) and (c) in Video S2 show the spatiotemporal evolution of V m for these cases. These animations, the pseudocolor plots of V m , the representative local time series of V m , the plots of the IBI, and the power spectra show that the initial conditions IC2 and IC3 lead, respectively, to spatiotemporal chaos and spiral turbulence (ST), with a single spiral meandering chaotically, and broken spirals, respectively, in the simulation domain.
Figures S9A-G, S10A-G, and S11A-F (Supplementary Material S1) show, respectively, the TNNP04 analogs of the TP06  Figures S9A-G, S10A-G, S11A-F (Supplementary Material S1) we conclude that the spatiotemporal evolution of V m in the TNNP04 model, without PD, is similar to, but not identically the same as, that in the TP06 model for the initial conditions IC1, IC2, and IC3. One difference is that, in the TNNP04 model, we have a Z-type, spiral-tip trajectory in Figures S10C,G (Supplementary Material S1), whereas, for the same initial condition, we have an open spiral-tip trajectory (Figures S7C,G in the Supplementary Material S1) in the TP06 model. This shows that spiral-wave dynamics in these two models, without PD, depends sensitively on the ionic details of these models and on the initial conditions.

SPIRAL WAVES WITH PD
We present systematic studies of spiral-wave dynamics here by using IC1, IC2, and IC3 initial configurations in the presence of PD. In the TP06 model, these initial configurations lead, respectively, to (a) an RS state with a roughly circular spiral-tip trajectory, (b) a single meandering spiral with turbulence (we refer to this as SMST henceforth), and (c) multiple-spiral turbulence (MST) with broken spiral waves in the absence of PD, as we have described above. For the TNNP04 model the analogs of these states are (a) RSC, a state with a rotating spiral whose tip trajectory is roughly circular, (b) RSZ, a state with a rotating spiral whose tip trajectory is roughly Z-type, and (c) an MST state.
We first consider the time evolution of IC1 for the TP06 model in the presence of PD, for which we deform the medium periodically along both x and y directions, with amplitudes and frequencies in the ranges 0.1 ≤ A x = A y ≤ 0.5 and 1.0 Hz ≤ f x = f y ≤ 7.0 Hz, respectively.
In Figures 4A-D we show pseudocolor plots of V m at time t = 4 s for (a) A x = A y = 0.1, f x = f y = 1.0 Hz, (b) A x = A y = 0.1, f x = f y = 3.0 Hz, (c) A x = A y = 0.1, f x = f y = 5.0 Hz, and (d) A x = A y = 0.1, f x = f y = 7.0 Hz, respectively. The RS state, which we obtain in the absence of PD, does not evolve into an MST state in cases (a), (b) and (c); however, in case (d) the spiral arm splits into multiple spirals to yield an MST state with mild spatiotemporal chaos, in so far as the dominant spiral does not break down but continues to evolve somewhat like a mother rotor (Samie and Jalife, 2001;Chen et al., 2003;Wu et al., 2004;Ideker and Rogers, 2006); the animations (a1), (b1), (c1), and (d1) in Video S3 show the spatiotemporal evolution of V m for these cases in the time interval 0 s ≤ t ≤ 4 s; this video uses 10 frames per second (fps) and each pseudocolor plot of V m is separated from its predecessor by 8 ms. The spiral-tip trajectories, which follow from this spatiotemporal evolution, are shown in Figures 3E-H for 3.6 s ≤ t ≤ 4 s (in Figure 4H we give the tip trajectory for the main, central spiral in Figures 4A-D); these tip trajectories are nearly circular with radii l c 18 mm, but, as we show below, the temporal evolution of V m is different in these cases. To examine this evolution, we obtain the local time series of V m (x, y, t), from the representative point (x = 125 mm, y = 125 mm) (the asterisks in Figures 4A-D), and therefrom the plots of the IBI shown in  Figure 4). From Figures 4I-L we see that the IBI displays a slight upward trend; this implies that, although the temporal evolution is nearly periodic, there is a gentle drift, toward lower frequencies, in the rotation rate of the dominant spiral. Furthermore, there are small oscillations in the IBI in Figure 4I (a 5-cycle), (j)(a 3-cycle), and (l)(a 2-cycle), but not in Figure 4K  Similarly, we study the dependence of spiral-wave dynamics on the amplitudes A x and A y of the PD, with the frequencies f x = f y held at a fixed value. In Figure 5 we show the pseudocolor plots of V m at time t = 4 s for (a) Hz. The spiral wave does not split into multiple spirals for these representative values of the amplitudes and frequencies. The animations (a1), (a2), (a3), and (a4) in Video S3 show the spatiotemporal evolution of these spiral waves for the interval 0 s ≤ t ≤ 4 s. To examine this evolution, we obtain the local time series of V m (x, y, t), from the representative point (x = 125 mm, y = 125 mm) (the asterisks in Figures 5A-D) and the corresponding tip trajectories of spiral waves, in the time interval 3.6 s ≤ t ≤ 4 s (blue lines with black points in Figures 5E,H The discrete peaks in E(ω) appear at the following frequencies: (M) ω 1 = 4.75 Hz, ω 2 = 9.5 Hz, ω 2 = 14.25 Hz, (N) ω 1 = 4.75 Hz, ω 2 = 9.5 Hz, ω 2 = 14.25 Hz, and small peaks at ω 1 = 3 Hz, ω 2 = 7.75 Hz, ω 3 = 11 Hz, ω 4 = 12.5 Hz, ω 5 = 15.75 Hz, (O) ω 1 = 4.75 Hz, ω 2 = 9.5 Hz, ω 2 = 14.25 Hz, and (P) ω 1 = 4.75 Hz, ω 2 = 9.5 Hz, ω 2 = 14.25 Hz. In (I-L) we see that the IBI shows a slight upward trend; this implies that, although the temporal evolution is nearly periodic, there is a slight drift, toward lower frequencies, in the rotation rate of the dominant spiral; note also the mild oscillations in the IBI in (I) a 5-cycle, (J) a 3-cycle, and (L) a 2-cycle, but not in (K) a 1-cycle; the natures of these oscillations and their cycle lengths are confirmed by the Poincaré-type return maps, shown in (Q-T), respectively; in these return maps, successive points are connected by lines. discrete peaks in E(ω) appear at the fundamental frequency ω f 4.75 Hz and at the frequencies listed in the caption of Figure 5; these peaks indicate that, in Figures 5M,N, we also have some high-order cycles; the broad-band power spectra in Figures 5O,P provide evidence for spiral turbulence with a meandering spiral (SMST). In Figures 5Q-T, we show Poincaré-type return maps which we obtain from the IBI plots in Figures 5I-L; in these maps successive points are connected by lines. These plots give  Figures 4A-D) and therefrom the plots of the IBI shown in (I-L), the power spectra E(ω) in (M-P), and the Poincaré-type return maps, which we obtain from the IBI plots and which show (Q) a 5−cycle, (R) a 5−cycle, (S) chaotic behavior, and (T) chaotic evolution; discrete peaks in E(ω) appear at the following frequencies: (M) ω 1 = 4.75 Hz, ω 2 = 9.5 Hz, ω 2 = 14.25 Hz, (N) ω 1 = 4.75 Hz, ω 2 = 9.5 Hz, ω 2 = 14.25 Hz and small peaks at ω 1 = 3.75 Hz, ω 2 = 8.5 Hz, ω 3 = 13.25 Hz, ω 4 = 15 Hz, ω 5 = 17.75 Hz, (O) ω 1 = 4.75 Hz, ω 2 = 9.5 Hz, and (P) ω 1 = 4.5 Hz, ω 2 = 9.5 Hz.
additional evidence for five cycles in Figures 5I,J We focus next on the types of ST states that we obtain, with PD applied along both x and y axes, when we start with the IC1 initial condition. In and (x = 50 mm, y = 50 mm), respectively. These pseudocolor plots and animations of V m and the plots of the IBI and power spectra show that we have, roughly speaking, three types of ST states with (a) multiple spirals (Figure 6A), (b) a stable spiral core with broken spiral arms (Figure 6B), and (c) a single, dominant, meandering spiral ( Figure 6C); the second case (b) displays a coexistence of a quasiperiodic and an ST state because of the dominant spiral at the center and the broken spirals generated from its arm. Such coexistence behaviors have been observed in both computational (Xie et al., 2001;Fenton et al., 2002;Cherry and Fenton, 2008) and experimental studies (Nash et al., 2006;Massé et al., 2007), which include in vivo experiments.
We also obtain quiescent (Q) states with no spirals because of the absorption of spiral waves at the boundaries as shown in the pseudocolor plots of V m in Figure S12 in the Supplementary Material S1 and in Video S3.
To illustrate the rich variety of spatiotemporal patterns, we summarize our results for the TP06 model, with the initial condition IC1, by presenting a selection of pseudocolor plots of V m in Figure S13(a1-d5) in the Supplementary Material S1 (for parameter sets see Table 1). The animations in Video S3 show the spatiotemporal evolution of V m for these cases in the time interval 0 s ≤ t ≤ 4 s. To examine this evolution, we obtain the local time series of V m (x, y, t), from the representative points (x = 125 mm, y = 125 mm); these are shown in Figure S14 (Supplementary Material S1); from these local time series, we obtain the plots of the IBI ( Figure S15 in the Supplementary Material S1) and the power spectra ( Figure S16 in the Supplementary Material S1).
The counterparts of Figures S13-S16 in the Supplementary Material S1, for initial conditions IC2 for IC3, are given, respectively, in Figures S17-S24 in the Supplementary Material S1.
For the initial conditions IC2 and IC3 the analogs of the animations in Video S3 are given, respectively, in Videos S4, S5. For IC2, with PD along both axes and different values of the amplitude and the frequency, we examine the time series of V m (x, y, t), from a representative point in the simulation domain ( Figure S18 in the Supplementary Material S1), the plots of the IBI ( Figure  S19 in the Supplementary Material S1) and the power spectrum E(ω) ( Figure S20 in the Supplementary Material S1), and the spatiotemporal evolution of V m (given by the animations in Video S4) and conclude therefrom that, in this case, we obtain either (a) a Q state with no spirals [see animations (a4), (a5), (b4), (b5), (d4) and (d5) in Video S4] or (b) an MST state with broken spiral waves (see the remaining animations in Video S4). A similar analysis, for IC3 and PD along both x and y axes, based on time series of V m ( Figure S22 in the Supplementary Material S1), plots of the IBI ( Figure S23 in the Supplementary Material S1), the power spectrum ( Figure S24  The TNNP04 model with PD also exhibits a rich variety of spatiotemporal patterns with spiral waves like the TP06 model as we discuss in the Supplementary Material S1 (see Figures  S25-S36). These figures and the associated Videos S6-S8 show that spiral-wave dynamics with PD in the TNNP04 model is quantitatively different from, but qualitatively similar to, that in the TP06 model. These differences arise because of the differences in the calcium-ion dynamics in these models; such dynamics can play an important role in the spatiotemporal evolution of spiral waves (Weiss et al., 2005;Ter Keurs and Boyden, 2007).
We have discussed spiral-wave dynamics in TP06 and TNNP04 models in the presence of PD along both x and y directions, with initial conditions of types IC1, IC2, and IC3. We have also carried out systematic simulations of spiral-wave dynamics in both these models, with PD along only one (say x) direction. Here too our results are, in the main, qualitatively similar to those we have presented above. Of course, there is anisotropic diffusion if the PD is only along one direction. However, Q, RS, and ST states appear; an overview of their spatiotemporal evolution is given in Figures  S37-S42 in the Supplementary Material S1.

SUPPRESSION OF SPIRAL WAVES
One of the principal goals of our extensive numerical studies of spiral-wave dynamics in the TP06 and TNNP04 models with PD is to understand its role in enhancing or suppressing spiralwave turbulence; this is an important step in developing an effective, low-amplitude suppression technique for the elimination of turbulence with single or multiple spirals. So far, various low-amplitude suppression algorithms have been developed to eliminate spiral waves in monodomain mathematical models of cardiac tissue; in these algorithms the control pulses are applied in several ways. These include the following: (a) periodic stimulation at a point (Zhang et al., 2003;Yuan et al., 2005); (b) a line stimulus that must be applied to one of the boundaries (Tang et al., 2008;Miguel et al., 2009); (c) an array of low-voltage control pulses, which must be swept over the simulation domain (Sinha and Sridhar, 2007;Sridhar and Sinha, 2008); or (d) the meshbased, low-amplitude suppression scheme we describe (Sinha et al., 2001;Pandit et al., 2002). We have provided an overview of such low-amplitude suppression schemes, in the absence of PD, in earlier studies (Shajahan et al., 2009;Nayak, 2013); the most successful of these is based on a mesh-based suppression algorithm; this suppression scheme (Shajahan et al., 2009;Majumder et al., 2011a;Nayak et al., 2013) can suppress spiral waves of electrical activation even in the presence of conduction, ionic, and fibroblast heterogeneities (Shajahan et al., 2009;Majumder et al., 2011a;Nayak et al., 2013). We now investigate the efficacy of this meshbased suppression scheme for both TP06 and TNNP04 models in the presence of PD.
In this mesh-based suppression scheme, we apply a current pulse of amplitude 75 pA/pF for 0.2 s over a mesh that divides our square simulation domain with L = 256 mm into 64 square cells of side l = 32 mm each; this pulse makes the links of the mesh refractory and, thereby, effectively imposes Neumann boundary conditions for any block inside the mesh; therefore, spiral waves inside a block are absorbed on the links of the mesh that bound the block. We have also extended this mesh-based scheme to one that uses control pulses on a set of parallel lines; in this linebased scheme, we apply a current pulse of amplitude 125 pA/pF for 0.6 s over a set of parallel lines separated from each other by l = 32 mm. As we show in the Supplementary Material S1 (see Figures S43, S44), both these schemes succeed in suppressing spiral-wave turbulence in the TP06 and TNNP04 models without PD; the line-based scheme uses a higher amplitude for the control pulse and a longer duration of application than the meshbased one because the former has fewer control-pulse segments than the latter. Furthermore, we show (see Figures S45-S48, in the Supplementary Material S1) that the line-based scheme works with PD only if the PD is applied along one spatial direction.
A minor modification of our line-based suppression scheme suppresses spiral-wave turbulence: we use a rectangular-meshbased control scheme, in which we add a few control lines perpendicular to the parallel lines of the line-based suppression scheme. We present a comparison of spiral-wave suppression by low-amplitude pulses on square, line, and rectangular suppression meshes in the TP06 model, with PD along both x and y directions: We impose PD along both x and y directions with the illustrative amplitudes A x = A y = 0.3 and frequencies f x = f y = 5 Hz for the initial configurations IC1, IC2, and IC3 (pseudocolor plots of V m in Figures 7A,E,I, respectively). We apply the following control pulses: (1) pulses with amplitude 75 pA/pF for t = 0.2 s over a square mesh (Figures 7B,F,J), with each square block of side l = 32 mm; (2) pulses with amplitude 125 pA/pF for t = 0.6 s over a line mesh (Figures 7C,G,K), with inter-line spacing l = 32 mm; (3) pulses with amplitude 125 pA/pF for t = 0.6 s over a rectangular mesh (Figures 7D,H,L), with block sides l x = 32 mm and l y = 64 mm. These pseudocolor plots of V m and the associated animations in Video S9 show that such spiral-wave states, with IC1, IC2, and IC3 initial conditions, are suppressed by both square-and rectangular-mesh suppression but not by linemesh suppression. Our rectangular-mesh suppression scheme is a FIGURE 7 | Comparison of spiral-wave suppression by low-amplitude pulses on square, line, and rectangular control meshes in the TP06 model, with PD along both x and y directions. We impose PD along both x and y directions with the illustrative amplitudes A x = A y = 0.3 and frequencies f x = f y = 5 Hz for the initial configurations IC1, IC2, and IC3 (pseudocolor plots of V m in (A,E,I), respectively). We apply the following control pulses: amplitude 75 pA/pF for t = 0.2 s over a square mesh (B,F,J), with each square block of side l = 32 mm; amplitude 125 pA/pF for t = 0.6 s over a line mesh (C,G,K), with inter-line spacing l = 32 mm; amplitude 125 pA/pF for t = 0.6 s over a rectangular mesh (D,H,L), with block sides l x = 32 mm and l y = 64 mm. These pseudocolor plots of V m and the associated animations in Video S9 show that these spiral states, with IC1, IC2, and IC3 initial conditions, are suppressed by both square-and rectangular-mesh control but not line-mesh control.
significant improvement over the square-mesh one because it uses fewer control lines than the latter. The results of similar studies for the TNNP04 model are given in Figure S49 (Supplementary Material S1).

DISCUSSION AND CONCLUSION
We have carried out detailed and systematic numerical studies of the effects of periodic deformation (PD) on spiral-wave dynamics in ionically realistic mathematical models for cardiac tissue by introducing PD in the recently developed TP06 and TNNP04 mathematical models for human ventricular tissue (Ten Tusscher et al., 2004;Ten Tusscher and Panfilov, 2006), in which we include deformation as in Zhang et al. (2004Zhang et al. ( , 2006 and Chen et al. (2008). We also investigate, in 2D simulations with PD, the efficacies of square-, rectangular-, and line-mesh-based, low-amplitude suppression techniques to eliminate spiral-wave turbulence in these models for cardiac tissue (Sinha et al., 2001;Pandit et al., 2002;Shajahan et al., 2009;Majumder et al., 2011a).
We have first considered simulations in cable-type domains, which are ideally suited for the calculation of θ and λ. We find that PD leads to a periodic, spatial modulation of θ and a temporally periodic modulation of λ; the degrees of these modulations depend on the amplitude A x and frequency f x of the PD. To the best of our knowledge, such modulations have not been quantified in any earlier study, although a few studies (Zhang et al., 2006;Chen et al., 2008) have suggested, in the context of spiral waves, that such modulations can arise because of a Doppler-type effect (Fenton et al., 2002).
We have considered three types of initial spiral-wave configurations, IC1, IC2, and IC3, for both the TP06 and TNNP04 models. In the TP06 model, these configurations evolve, respectively, to (a) an RSC state, (b) an SMST state, and (c) an MST state, in the absence of PD; in the TNNP04 model they evolve, respectively, to (a) an RSC state, (b) an RSZ state, and (c) an MST state, in the absence of PD. We have used such initial conditions because various experimental and computational studies (Damle et al., 1992;Davidenko et al., 1992;Pertsov et al., 1993;Gray et al., 1995;Gray and Jalife, 1996;Beaumont et al., 1998;Gray et al., 1998;Witkowski et al., 1998;Qu et al., 2000a;Fenton et al., 2002;Ten Tusscher and Panfilov, 2006;Shajahan et al., 2009) have suggested that spiral-wave dynamics in cardiac tissue can lead to (a) a stable rotor (Davidenko et al., 1992;Pertsov et al., 1993;Beaumont et al., 1998), as in our RSC state, (b) a single, meandering rotor whose time series is chaotic (Gray et al., 1995;Gray and Jalife, 1996;Qu et al., 2000a), as in our SMST state, and (c) multiple rotors, which yield a state with spatiotemporal chaos (Damle et al., 1992;Gray et al., 1998;Witkowski et al., 1998;Qu et al., 2000a), as in our MST state. Thus, our initial conditions, IC1, IC2, and IC3, lead to the three major types of spiral-wave evolutions, and slight variants thereof (e.g., RSZ), which have been seen in earlier studies and whose evolution we study now with PD. This shows that spiral-wave dynamics in these two models, without PD, depends sensitively on the ionic details of these models and on the initial conditions. A rich variety of spiralwave behaviors result when we add PD to the TP06 and TNNP04 models, which are quantitatively different from, but qualitatively similar to, that in the TP06 model. These differences arise because of the differences in the calcium-ion dynamics in these models; such dynamics can play an important role in the spatiotemporal evolution of spiral waves (Weiss et al., 2005;Ter Keurs and Boyden, 2007). Our principal findings here can be summarized as follows: In the presence of PD, an RS state may show (a) periodic behavior with high-order cycles in time, (b) temporally quasiperiodic (QP) evolution, (c) a state with spiral-wave turbulence, or (d) a quiescent state Q. For an ST state, which can be of SMST or MST types, PD can either leave the system in an ST state or make it evolve to a Q state, in which all spirals either annihilate each other or are absorbed at the boundaries of the simulation domain. Precisely which one of these states is obtained depends sensitively on our initial conditions and on the PD parameters A x , A y , f x , and f y . This sensitive dependence on parameters is a hallmark of extended dynamical systems that show spatiotemporal chaos (Shajahan et al., 2007(Shajahan et al., , 2009). Thus, our study systematizes the effects of PD on spiral-wave dynamics and turbulence in two, biophysically realistic mathematical models for cardiac tissue; and it complements earlier studies of spiral-wave dynamics, in such models, that have concentrated on the dependence of such dynamics on ion-channel and electrophysiological properties (Karma, 1994;Qu et al., 1999Qu et al., , 2000a and on conduction (Xie et al., 1998;Ten Tusscher and Panfilov, 2003;Shajahan et al., 2007Shajahan et al., , 2009) and ionic inhomogeneities (Shajahan et al., 2009;Majumder et al., 2011a,b). By using the biophysically realistic TP06 and TNNP04 models for cardiac tissue, our study generalizes the work of Zhang et al. (2004Zhang et al. ( , 2006 and Chen et al. (2008) on spiral-wave instabilities in a simple, two-variable model for cardiac tissue, which is subject to PD.
Moreover, our studies have used three types of spiral-wave initial configurations to examine, via extensive and systematic numerical calculations, the transitions between different states of our system as the amplitude and frequency of the PD are varied. Our work extends significantly earlier studies of PD (Zhang et al., 2004(Zhang et al., , 2006Chen et al., 2008) and mechanical deformation (Panfilov et al., 2005(Panfilov et al., , 2007Weise et al., 2011). Zhang et al. (2004) have studied the instability of a spiral wave of electrical activation by introducing, in a simple, two-variable, FitzHugh-Nagumo-type model (Aliev and Panfilov, 1996) for cardiac tissue, the possibility of periodic, temporal oscillations in the diffusion constant. Their study shows that the resulting deformation can lead to a transition from a stable, RS state to an ST state with multiple spirals. In another study Zhang et al. (2006) have shown that such an ST state can be driven into a quiescent state with no spirals, if the oscillation frequency of the PD is chosen to be close to the characteristic frequency of the spiral wave in the RS state of the system. Chen et al. (2008) have studied the effects of PD in the two-variable, Bär model (Bär and Eiswirth, 1993) on spiral-wave dynamics by varying the parameter , which sets the time scale of the slow variable in this model, and the amplitude A and frequency f of the PD; their study shows that the RS-ST transition can be effected by changing , A, and f suitably. Panfilov et al. (2007) have shown that mechanical deformation can either (a) induce a spiral wave to drift or (b) break up spiral waves and thus lead to complex spatiotemporal patterns in the three-variable Fenton-Karma model (Fenton and Karma, 1998) for cardiac tissue; the model that (Panfilov et al., 2007) use for deformation is different from, and more realistic than, the one used in Zhang et al. (2004Zhang et al. ( , 2006 and Chen et al. (2008) in so far as it includes a stretch-activated current, which accounts for the mechano-electrical feedback in cardiac tissue, whose stress tensor controls the deformation; their study shows that rotating spirals become unstable both because of the stretch-activated current and the deformation of the tissue. In a related study, which also includes stretch-activated currents, Weise et al. (2011) have shown that such deformation can lead to pacemaker activity in a discrete version of the two-variable, Aliev-Panfilov, reactiondiffusion model (Aliev and Panfilov, 1996). Note that the studies in Zhang et al. (2004Zhang et al. ( , 2006, Panfilov et al. (2007), Chen et al. (2008), and Weise et al. (2011) have used only a particular type of spiral-wave configuration in their two-variable models, with either PD or mechanical deformation. Moreover, Zhang et al. (2004Zhang et al. ( , 2006, Panfilov et al. (2005Panfilov et al. ( , 2007, Chen et al. (2008), and Weise et al. (2011) have focused on simple, two-variable, mathematical models for cardiac tissue. Thus, these studies cannot address spiral-wave dynamics in such tissue at the detailed ionic level we consider in our work by using state-of-the-art, ionically realistic mathematical models for ventricular tissue. Furthermore, the authors of Zhang et al. (2004Zhang et al. ( , 2006 and Chen et al. (2008) study the effects of PD on spiral-wave dynamics for a limited set of initial conditions. For example, in Zhang et al. (2004), the authors have studied the behavior of a single rotating spiral (RS) in the presence of PD; the authors of Zhang et al. (2006) have used a broken-spiral state as an initial configuration to study the elimination of spirals from the system in the presence of PD; in Chen et al. (2008), the authors have studied an RS initial configuration and its spatiotemporal evolution with PD. The authors of Panfilov et al. (2005Panfilov et al. ( , 2007 and Weise et al. (2011) have used an RS initial configuration to examine the effect of mechanical deformation on this RS state; they have not investigated the transitions between different spiral-wave states. None of these studies have carried out the detailed numerical investigations of spiral-wave dynamics that we present in our work, which considers a variety of initial conditions. Furthermore, we have shown that square-and line-meshbased, low-amplitude suppression schemes eliminate spiral-wave turbulence in both the TP06 and TNNP04 models in the absence of PD; this line-based suppression scheme is a significant improvement over the square-mesh suppression scheme of Sinha et al. (2001), Pandit et al. (2002), Shajahan et al. (2009), andMajumder et al. (2011a) because it has fewer control lines. However, we have found that the line-based scheme works with PD only if the PD is applied along one spatial direction. We have then shown that a minor modification of our line-based suppression scheme can suppress spiral-wave turbulence: in particular, we introduce a rectangular-mesh-based suppression scheme, in which we add a few control lines perpendicular to the parallel lines of the line-based suppression scheme; this rectangular-mesh scheme is also a significant improvement because it uses fewer control lines than the one based on a square mesh. We would like to emphasize here that no spiral-wave suppression scheme has been tried in the presence of PD hitherto; our study is the first to address this issue.
The formation of patterns in reaction-diffusion type system with various types of flows have been investigated in Biktashev et al. (1998), Leconte et al. (2003), Kuptsov et al. (2005), and Ermakova et al. (2009); in particular, break up of spiral excitation waves has been observed in a moving excitable medium as suggested in Biktashev et al. (1998); these studies have shown that linear shear flow can cause spiral-wave breakup in an excitable medium. Recent studies in Yoshida (2010) and Yashin et al. (2012a,b) have investigated pattern formation in a gel medium that can be oscillated mechanically; the effect of these mechanical oscillations on the underlying spatiotemporal chemical oscillations, because of a Belousov-Zhabotinsky (BZ) reaction, can be studied in such systems. Our work here presents the cardiac-tissue analogs of such chemical-oscillation studies.
Our study has used methods of non-linear-dynamics to evaluate the effects of periodic deformation on waves of excitation, both plane waves and spiral waves, and has obtained therefrom important physiological implications. In particular, we have shown, for the first time, how spiral-wave turbulence can be suppressed even in the presence of such deformation. This has potential applications in defibrillation because spiral-wave turbulence in our mathematical models is the analog of ventricular fibrillation (Davidenko et al., 1992;Pertsov et al., 1993;Gray et al., 1998;Jalife et al., 1998). Furthermore, our study has immediate implications for experiments on cell cultures with cardiac cells.
Various stretching devices have been developed to control the contraction and expansion of a cell (Barbee et al., 1994;Huang et al., 2010;Wang et al., 2010) and a layer of cells in culture Sotoudeh et al., 1998;Waters et al., 2001;Winter et al., 2002;Rana et al., 2008). In these devices, both uniaxial and biaxial (Pang, 2009) stretching methods can be used to deform substrates; moreover, in some of these devices, the stretching can be applied in a cyclic manner at different frequencies. This stretching-induced deformation of the substrate leads, in turn, to the deformation of a layer of culture cells that are attached to the substrate. Examples of such studies include the following: The study in Barbee et al. (1994) has measured the strain that develops in cultured vascular smooth muscle cells when they are deformed by the stretching of a substrate to which they adhere. The authors of Lee et al. (1996) have used a device, which applies homogeneous, equibiaxial strains of 0-10% to a cell-culture substrate, to verify quantitatively the transmission of substrate deformation to a 2D sheet of cultured cardiac cells. The studies in Sotoudeh et al. (1998) have used endothelial cells in tissue culture, on a silicon elastic membrane, and have designed an apparatus that allows for the control of the magnitude and frequency of the dynamical stretching that is applied uniformly to these cells to produce equibiaxial dynamical stretches, with area changes ranging from 0% to 55% and frequencies ranging from 0 to 2 Hz. The authors of Waters et al. (2001) have developed a system for the imposition of cyclic biaxial strain to stretch cultured pulmonary epithelial cells; similar techniques have been used in Winter et al. (2002) to study the effects of strain in cell cultures and in vitro experiments. The authors of Rana et al. (2008) have studied the response to such stretching in cultured neonatal rat atrial cardiomyocytes by using a device that can impose homogeneous equibiaxial deformation. Other recent studies include those of Huang et al. (2010) and Wang et al. (2010), which have studied the mechanical activities of living cell, fiber, and tissue by applying both equiaxial and uniaxial deformation, and recording the dynamics of the response of these systems by using highresolution imaging techniques; the former experiment has used fibroblasts and the latter endothelial cells in culture. We suggest that such experimental studies of the responses of cell cultures to an applied stress can be easily generalized to study the types of problems we have concentrated on here. In particular, by imposing a periodic deformation on cardiac tissue or cell cultures, experiments should be able to verify the predictions we have made, on the basis of our in silico studies, about the modulations of θ and λ in the presence of PD and the effects of PD on spiralwave dynamics, which we have discussed in detail in the previous Section.
We end with some limitations of our model. The first limitation is that our model does not include stretch-activated currents (Panfilov et al., 2005(Panfilov et al., , 2007Weise et al., 2011) explicitly; but V m depends on PD and all ionic currents depend on V m , so PD affects all such currents implicitly. Next, the PD in our model affects the electrical activation of our medium but it is not, in turn, affected by this activation; by contrast, the model for mechanical deformation used in Panfilov et al. (2005Panfilov et al. ( , 2007 and Weise et al. (2011) allows for electrical feedback to affect such deformation; our model does not include soft-tissue mechanics, which can be incorporated in mathematical models for cardiac tissue by including stress and strain tensors, from elasticity theory (Nash and Hunter, 2000;Nash and Panfilov, 2004;Keldermann et al., 2009), as in the studies of Panfilov et al. (2005Panfilov et al. ( , 2007 and Weise et al. (2011); however, these studies use only a twovariable model for cardiac tissue and not the ionically realistic TP06 or TNNP04 models that we employ. Moreover, because of the absence of detailed ion-channel dynamics, the simple, twovariable models for cardiac tissue, which have been used in the studies of Panfilov et al. (2005Panfilov et al. ( , 2007 and Weise et al. (2011), do not account for the effects of deformation on ion-channel activity and the intracellular calcium concentration as suggested in Cherubini et al. (2008), Pathmanathan andWhiteley (2009), andAmbrosi et al. (2011). In spite of the simplicity of our model for PD, our study captures various features of spiral-wave dynamics that have been observed in models that include stretch-activated currents (Panfilov et al., 2007); in particular, our model displays spiral-wave breakup because of PD. Thus, this result of ours is robust. The only qualitative effect that our study misses is deformation-induced pacemaker activity, for which it has been argued (Kohl et al., 1999;Panfilov et al., 2005;Weise et al., 2011) that stretch-activated currents are essential. To the best of our knowledge, our elucidation of the effects of PD on spiral-wave dynamics in mathematical models for cardiac tissue, although simple in its modeling of PD, is the first study that explores the effects of PD on spiral-wave dynamics in ionically realistic mathematical models for ventricular tissue. A complete study of a realistic model for deformation, with stretch-activated currents, and such ionically realistic mathematical models lies beyond the scope of the present paper. While this paper was being prepared for publication, a new study appeared on a human-ventricular mathematical model with mechanical deformation (Weise and Panfilov, 2013). The authors of this study have used the same type of deformation that they have employed in their previous investigations of mechanical deformation (Weise et al., 2011); they have studied spiral-wave dynamics in the context of the drifting of a spiral wave as a function of G s (here, G s is the maximum conductance of the stretch-activated current). In their numerical studies they have found that, in a constantly stretched medium, there is an increase of the core size and period of a spiral wave, but no change in its rotational dynamics; in contrast, in the dynamically stretched medium, they observe spiral drift. We have found such behaviors, to some extent, in our periodically deformed medium; e.g., the initial condition IC1, which leads to a stable rotating spiral with a circular tip trajectory, in the absence of PD, can lead to a variety of behaviors, which can include slow spiral drift in the presence of PD. Our model does not include mechano-electrical feedback in a realistic way, as we have described above. Therefore, we have not attempted to study how different mechanical stimuli, other than the PD we consider, initiate or affect spiral-waves in our model; studies of other mechanical stimuli lie beyond the scope of our paper. In Weise et al. (2011) it has been noted that both electrical and mechanical stimuli can cause the formation of a pacemaker in cardiac tissue; and mechanical stimuli can translate the mechanical energy into an electrical stimulus, as argued in Janse et al. (2003) and Cooper et al. (2006). We use a monodomain description for cardiac tissue; and we do not use an anatomically realistic simulation domain (Panfilov and Keener, 1995;Trayanova and Tice, 2009), muscle-fiber orientation, and transmural heterogeneity (Majumder et al., 2011b(Majumder et al., , 2012; the inclusion of these features lies beyond the scope of this study. We note, however, that recent studies (Potse et al., 2006) have compared potentials resulting from normal depolarization and repolarization in a bidomain model with those of a monodomain model; these studies have shown that the differences between results obtained from a monodomain model and those obtained from a bidomain model are extremely small.

ACKNOWLEDGMENTS
We thank the Department of Science and Technology (DST), India, the University Grants Commission (UGC), India, Council for Scientific and Industrial Research (CSIR), India, Jawaharlal Nehru Centre for Advanced Scientific Research (JNCASR), India, and the Robert Bosch Centre for Cyber Physical Systems (IISc) for support, and the Supercomputer Education and Research Centre (SERC, IISc) for computational resources.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/fphys.2014. 00207/abstract Video S1 | Spatiotemporal evolution of plane waves in cable-type domains, with PD along the axial-direction of the cables, for the TP06 model, shown via pseudocolor plots of the transmembrane potential V m ; the time evolution is shown for 0 s ≤ t ≤ 2 s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S2 | Spiral-wave dynamics in the TP06 and TNNP04 models in the absence of PD; the time evolution of pseudocolor plots of the transmembrane potential V m is shown for 0 s ≤ t ≤ 4 s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S3 | Spiral-wave dynamics for the TP06 model in the presence of PD along both x and y directions with an initial condition of type IC1; the time evolution of pseudocolor plots of the transmembrane potential V m is shown for 0 s ≤ t ≤ 4 s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S4 | Spiral-wave dynamics for the TP06 model in the presence of PD along both x and y directions with an initial condition of type IC2; the time evolution of pseudocolor plots of the transmembrane potential V m is shown for 0 s ≤ t ≤ 4 s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S5 | Spiral-wave dynamics for the TP06 model in the presence of PD along both x and y directions with an initial condition of type IC3; the time evolution of pseudocolor plots of the transmembrane potential V m is shown for 0 s ≤ t ≤ 4 s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S6 | Spiral-wave dynamics for the TNNP04 model in the presence of PD along both x and y directions with an initial condition of type IC1; the time evolution of pseudocolor plots of the transmembrane potential V m is shown for 0 s ≤ t ≤ 4, s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S7 | Spiral-wave dynamics for the TNNP04 model in the presence of PD along both x and y directions with an initial condition of type IC2; the time evolution of pseudocolor plots of the transmembrane potential V m is shown for 0 s ≤ t ≤ 4 s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S8 | Spiral-wave dynamics for the TNNP04 model in the presence of PD along both x and y directions with an initial condition of type IC3; the time evolution of pseudocolor plots of the transmembrane potential V m is shown for 0 s ≤ t ≤ 4 s; we use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.
Video S9 | Comparison of spiral-wave suppression by low-amplitude pulses on square, line, and rectangular control meshes in the TP06 model, with PD along both x and y directions, for the time interval 0 s ≤ t ≤ 1 s; we depict pseudocolor plots of V m and use 10 frames per second (fps); in real time each frame is separated from the succeeding frame by 8 ms.