Scroll-Wave Dynamics in Human Cardiac Tissue: Lessons from a Mathematical Model with Inhomogeneities and Fiber Architecture

Cardiac arrhythmias, such as ventricular tachycardia (VT) and ventricular fibrillation (VF), are among the leading causes of death in the industrialized world. These are associated with the formation of spiral and scroll waves of electrical activation in cardiac tissue; single spiral and scroll waves are believed to be associated with VT whereas their turbulent analogs are associated with VF. Thus, the study of these waves is an important biophysical problem. We present a systematic study of the combined effects of muscle-fiber rotation and inhomogeneities on scroll-wave dynamics in the TNNP (ten Tusscher Noble Noble Panfilov) model for human cardiac tissue. In particular, we use the three-dimensional TNNP model with fiber rotation and consider both conduction and ionic inhomogeneities. We find that, in addition to displaying a sensitive dependence on the positions, sizes, and types of inhomogeneities, scroll-wave dynamics also depends delicately upon the degree of fiber rotation. We find that the tendency of scroll waves to anchor to cylindrical conduction inhomogeneities increases with the radius of the inhomogeneity. Furthermore, the filament of the scroll wave can exhibit drift or meandering, transmural bending, twisting, and break-up. If the scroll-wave filament exhibits weak meandering, then there is a fine balance between the anchoring of this wave at the inhomogeneity and a disruption of wave-pinning by fiber rotation. If this filament displays strong meandering, then again the anchoring is suppressed by fiber rotation; also, the scroll wave can be eliminated from most of the layers only to be regenerated by a seed wave. Ionic inhomogeneities can also lead to an anchoring of the scroll wave; scroll waves can now enter the region inside an ionic inhomogeneity and can display a coexistence of spatiotemporal chaos and quasi-periodic behavior in different parts of the simulation domain. We discuss the experimental implications of our study.


Introduction
Cardiac arrhythmias, such as ventricular tachycardia (VT) and ventricular fibrillation (VF), are among the leading causes of death in the industrialized world. There is growing consensus [1][2][3][4][5][6][7] that such arrhythmias are associated with the formation of spiral and scroll waves of electrical activation in mammalian cardiac tissue; single spiral and scroll waves are believed to be associated with VT whereas their broken, turbulent analogs are thought to be the underlying cause of VF. Thus, the study of these waves and their eventual elimination from cardiac tissue is a problem of central importance in biomedical and biophysical science. The existence of spiral and scroll waves in mammalian cardiac tissue has been confirmed by a variety of in vitro and in vivo studies as reported, e.g., in Refs. [8][9][10][11][12][13][14][15][16][17]; and these have been complemented by in silico investigations of spiral-and scroll-wave dynamics in mathematical models of cardiac tissue [1][2][3][18][19][20][21][22][23][24]. In silico studies have a significant advantage over experimental ones when we consider rotating scroll waves or scroll-wave turbulence in three-dimensional cardiac tissue because we can visualize scroll-wave dynamics in three-dimensional computational domains in much more detail than is possible experimentally. Furthermore, we can tune parameters, such as the Calcium-channel conductance, the degree of muscle-fiber rotation, or the position or the size of a conduction inhomogeneity, with ease, and thus explore their effects on scroll waves. A detailed understanding of such effects is of paramount importance in controlling scroll-wave turbulence by low-amplitude electrical stimuli; the development of such control methods is the computational analog of low-amplitude defibrillation.
Strictly speaking a spiral wave is broken at its tip, and, in 3D, a scroll wave at its filament and the association of a single spiral wave with VT and multiple spiral waves with VF, is one possible scenario for these arrhythmias; VT can also exist because of ectopic activity and VF can arise from a chaotically drifting spiral wave, but without any further breakup.
Inhomogeneities in cardiac tissue affect spiral and scroll waves in many different ways. Experimental studies [25][26][27] have shown that conduction inhomogeneities can either lead to the anchoring or pinning of such waves or they can, occasionally, eliminate these waves completely [28]. The larger the obstacle the higher is the tendency of anchoring [25][26][27]; but, in some cases, wave pinning might not occur even if the obstacle is large; and such waves can get attached to small conduction inhomogeneities, but only infrequently [27].
Numerical simulations have also been used to investigate the effects of a variety of inhomogeneities on spiral-wave dynamics in mathematical models for cardiac tissue: For example, Ref. [29] contains an investigation of spiral waves in the two-dimensional (2D) Luo-Rudy Phase I (LRI) model in an annular domain, with parameters and initial conditions such that, in the absence of the inner hole of the annulus, spiral waves break up. As the radius of the hole is reduced, from 9:2 cm to 0, the simulation domain goes from a 1D ring to a 2D annulus, and, finally, to a homogeneous 2D domain: When the system is a 1D ring the wave goes around it; in the annular case the wave is a spiral pinned to the inner hole, if its radius is sufficiently large; as this radius decreases transitions occur first to a quasi-periodically rotating spiral wave and, eventually, to spiral-wave break up and spatiotemporal chaos [29] with the spirals not attached to the hole.
Reference [30] has studied the effects of a random distribution of non-excitable cells spiral waves in the Panfilov model; as the percentage of such non-excitable cells increases spiral-wave break up is suppressed.
Analytical and numerical studies of the interplay of spiral waves with a piecewise-linear obstacle appear in Ref. [31]. The principal results of this are the following: if the medium has a high excitability, a wave can go around the boundary of the inhomogeneity and continue thereafter as if no obstacle had been encountered; in the low-excitability case, the two free ends survive because they cannot rejoin beyond the obstacle so they develop into spiral waves that can break up subsequently; and, in addition to the excitability and wave front's local curvature, the shape of the inhomogeneity also affects the attachment of spiral waves to it.
Earlier studies from our group [2,3] have presented numerical investigations of spiral-wave dynamics in the presence of conduction and ionic inhomogeneities in the Panfilov, LRI, reduced-Priebe-Beuckelmann, and TNNP models in two-dimensional simulation domains; preliminary studies of scroll-wave interactions with such inhomogeneities have also been considered in three-dimensional simulation domains for the simple Panfilov model without fiber rotation. The principal qualitative result of these studies is that the dynamics of spiral waves depends very sensitively on the position, shape, size, and type (i.e., conduction or ionic) of inhomogeneity; for details we refer the reader to Refs. [2,3]; ionic inhomogeneities also lead to anchored spiral waves, but with dynamics that is considerably richer than in the case of conduction obstacles; e.g., we have found rotating spiral waves with period-4 and period-5 cycles and also states in which there is a coexistence of a periodically rotating spiral and with broken spiral waves, in the region of the ionic inhomogeneity. Our studies of the effects of inhomogeneities on scroll waves in mathematical models for cardiac tissue had not included fiber rotation in 3D simulation domains. This is the problem we address here.
Studies by several groups [20,23,[32][33][34][35][36] have shown that fiber rotation plays an important role in the spatiotemporal evolution of scroll waves in mathematical models of cardiac tissue. We summarize the results of these studies before we embark on our study of scroll waves with both fiber rotation and inhomogeneities. These studies have investigated the effects of fiber rotation on scroll-wave dynamics in simple two-variable, FitzHugh-Nagumotype models and also in ionic models such as the LRI and TNNP models.
The intramural rotation of cardiac muscle fibers makes wave propagation in ventricular tissue very anisotropic as noted, e.g., in Ref. [32] for simple FitzHugh-Nagumo-type models and for the Beeler-Reuter model. This anisotropy can twist scroll-wave filaments and thus modify their dynamics; if this twist is very large it can destabilize a single transmural filament, create vortex rings, which can, in turn, expand and create new filaments that move chaotically. These fibers also lead to a propagation anisotropy [20]: the conduction velocity along the fiber axis is approximately three times greater than it is across it. This anisotropy can lead to a twisting of scroll-wave filaments and break up of scroll waves even if the excitability is not low [37][38][39]. In Ref. [36] the effects of fiber rotation on scroll-wave dynamics have been studied numerically for the LRI model but without any inhomogeneities; with the original LRI-model parameters, these studies find that a minimum thickness is required for scroll-wave break-up, which proceeds via filament bending. The effects of the transmural rotation of fibers on scroll waves have also been studied for the LRI model in cubical domains [40] of thickness 2:9 mm particularly with a view to elucidating the effects of different fiberrotation rates on the dynamics of these waves; as this rate increases the susceptibility of the medium to re-entry decreases but complex spatiotemporal patterns can occur. The work of Ref. [35] elucidates the differences between the dynamics of scroll waves in cuboid and anatomically realistic simulation domains; this study includes fiber rotation, uses the TNNP model, but does not include inhomogeneities of the type we study here.
Our earlier work [2,3] has elucidated the effects of different types of inhomogeneities -conduction and ionic -on spiral-and scroll-wave dynamics in mathematical models of cardiac tissue but without fiber rotation. The studies we have summarized above have studied the effects of fiber rotation on scroll-wave dynamics on a variety of mathematical models for cardiac tissue but without conduction or ionic inhomogeneities. To the best of our knowledge, the study we present here is the first one to investigate the combined effects of fiber rotation and inhomogeneities on scroll-wave dynamics in the three-dimensional, detailed ionic TNNP model (principally with epicardial parameters).
The effect of transmural heterogeneities has been studied in detail in Refs. [41][42][43], in two-dimensional simulations, and in Ref. [35] in isotropic, orthotropic and anisotropic three dimensional geometries. However, such studies do not include the role of inhomogeneities, which have a profound effect on the dynamical behavior of spiral and scroll waves. In this paper we principally investigate the effects of such inhomogeneities in the epicardium, and, in three representative cases, in the whole thickness of the ventricular wall.
Our principal goal is to study the effects of cardiac fiber rotation (FR) on the dynamics of scroll waves in three-dimensional (3D) simulation domains, both without and with inhomogeneities, in the mathematical model of cardiac tissue due to ten Tusscher, Noble, Noble, and Panfilov (TNNP) [44]. By FR we mean the following: The fibers in each plane are rotated, with respect to those in adjacent planes, by an angle a; thus, the orientation of the fibers in a plane is specified by an angle h(z) that varies continuously with the transmural coordinate z. In a normal human heart the average total FR, from the endocardium to the epicardium, is^110 0 {120 0 . The distribution of this angle over the three-layered heart wall varies across mammalian species and, within a specie, from individual to individual. Our study here extends considerably our earlier work [3] on spiral-wave dynamics in the two-dimensional (2D) TNNP model; the principal qualitative result of this study was that such dynamics depends delicately on the positions, shapes, and sizes of conduction and ionic inhomogeneities.
In the 3D TNNP model that we study here with fiber rotation and with conduction or ionic inhomogeneities, we find the following qualitatively interesting results: 1. In addition to displaying a sensitive dependence on the positions, sizes, and types of inhomogeneities, scroll-wave dynamics also depends upon the degree of fiber rotation.
2. Scroll waves do not anchor easily to cylindrical conduction inhomogeneities with small radii; this anchoring tendency increases with the radius; this result is in agreement with the experimental findings of Refs. [27,[45][46][47][48][49]. We conjecture that our result is a three-dimensional manifestation of the findings of Ref. [50] in which the authors suggest that, if the radius of the closed trajectory of the tip of a spiral wave is less than the radius of an obstacle, then this wave is pinned to the obstacle; if not, then the wave can be removed by suitable pacing of the medium.
3. If the scroll-wave filament exhibits weak meandering, then there is a fine balance between the anchoring of this wave at the inhomogeneity and a disruption of wave-pinning by fiber rotation; by contrast, if this filament displays strong meandering, then anchoring is still suppressed by fiber rotation, but, in addition, the scroll wave can be eliminated from most of the layers only to be regenerated by a seed wave.
4. We find, furthermore, that ionic inhomogeneities, which we model by changing the values of certain conductances in a cylindrical region in our simulation domain, can also lead to an anchoring of the scroll wave; scroll waves can now enter the region inside an ionic inhomogeneity and the system can display a coexistence of spatiotemporal chaos and quasi-periodic behavior in different parts of the simulation domain.
The remaining part of this paper is organized as follows. The section entitled Results contains a detailed description of our results. It has two subsections: the first of these deals with our studies of the 3D TNNP model in the absence of inhomogeneities; the second subsection contains our studies of the 3D TNNP model with inhomogeneities. The section entitled Discussions contains a discussion of our results and the section labeled Methods is devoted to a description of the TNNP model and the numerical details of our simulations.

Results
We consider first scroll-wave dynamics in a rectangularparallelepiped simulation domain without FR or inhomogeneities in the 3D TNNP model. We begin with a simple scroll wave generated by stacking a set of spiral waves that we obtain from simulations on a two-dimensional (2D) TNNP model [3]. This simple initial scroll wave continues to rotate without significant modification of its original shape as can be seen from the volume rendering of V (x,y,z), at the representative time t~1:2s, in figure 1(a); the spatiotemporal evolution of V (x,y,z,t) for this case is shown in Video S1 in the supplementary material. Figure 1(b) shows a space-time (z,t) pseudocolor plot of V (x,y,z,t) at the representative point (x~45mm,y~45mm). Figure 1(c) shows the power spectrum E(v) of the time series V (x,y,z,t) from the representative point (x~45mm,y~45mm,z~1mm); from this we see that the scroll wave rotates quasi-periodically in the medium because the principal peaks in the power spectrum E(v) can be indexed as n 1 v 1f zn 2 v 2f , where n 1 and (x~45mm,y~45mm,z~1mm) are integers and v 1f^3 :5Hz and v 2f^7 :5Hz are two incommensurate frequencies.
We now investigate how the dynamical behavior of this scrollwave changes as we introduce FR and inhomogeneities in the medium. Subsection No inhomogeneities deals with the investigation of scroll-wave dynamics in the absence of localized inhomogeneities. Subsection Inhomogeneities deals with the combined effects of FR and inhomogeneities on scroll-wave dynamics; this subsection has two parts.

No inhomogeneities
Figures 2(a), (b), (c), and (d) show, at the representative time t~1:2s, volume renderings of V (x,y,z) for the FR angles Dh~0 0 ,10 0 ,30 0 and 60 0 , respectively; the spatiotemporal evolution of V(x,y,z,t) for the representative case of Dh~30 0 is shown in Video S2 in the supplementary material; the videos for the cases Dh~0 0 ,10 0 ,30 0 and 60 0 , are available at http://www.physics.iisc. ernet.in/*rahul/new_movies.html. Figures 2(e), (f), (g), and (h) show space-time (z,t) pseudocolor plots of V (x,y,z,t), at the representative point (x~45mm,y~45mm) for the same values of Dh; and figures 2(i), (j), (k), and (l) show the power spectra E(v) of the time series V (x,y,z,t) from the representative point (x~45 mm,y~45 mm,z~2mm). By comparing figure 1(a) with figures 2(a)-(d) (or Video S1 with the videos for the cases with Dh~0 0 ,10 0 ,30 0 and 60 0 ) we see that the most dramatic effect of FR is a distortion of the scroll-wave relative to its form in the completely homogeneous case; this distortion arises from the difference between D 11 and D 22 in the diffusion tensor (see the section labeled Methods); however, for the case with no fiber rotation, the TNNP model, with the tensor D, can be mapped onto the homogeneous case by suitable rescaling of the linear dimensions along and perpendicular to the fiber orientation.
A comparison of the space-time plots in figures 2(e) and (f) and of the power spectra in figures 2(i) and (j) shows that the spatiotemporal behaviors of the scroll waves for Dh~0 0 and Dh~10 0 are very similar. Deviations from the behavior for Dh~0 0 become noticeable for Dh~30 0 and very significant for Dh~60 0 ; the rotation frequencies remain approximately the same [see figures 2 (i)-(l)]; but the bowing in the space-time plots of figures 2(g) and (h) indicate that, as Dh increases, the spiral sections of the scroll wave are such that the one in the central plane z~10 (1 mm) leads those in the planes above and below it. This bowing is a manifestation of the transmural bending of the filamentary core of the scroll wave; and this bending can be seen clearly in Video S2, by looking at the curvature of the patterns on the boundaries that are normal to the (x,y) plane, and in the V~{30 mV isosurface video, Video S3 (see supplementary material) in which we have zoomed in to the region near the core of the scroll wave. Such a bending of the scroll-wave filament has been seen earlier in a numerical study of the Luo-Rudy I (LRI) model [36] that includes FR. In the LRI model the bending of the filament increases with Dh, as in the TNNP model we study; but, for sufficiently large Dh, the filament of the LRI scroll wave breaks. We have not found such filament breaking in the TNNP model, with the standard parameters [44] and for the range of values of Dh that we have used, even when we have increased the thickness of the simulation domain from z~2 mm to z~9 mm. For We find, however, that the 3D TNNP model does show scrollwave break-up if we change some parameters. In particular, if (a) we reduce the L-type Ca 2z conductance G CaL to one-fourth of its maximum value or (b) we increase the plateau Ca 2z conductance G pCa by a factor of 5, we obtain such break-up. In Figure 4 we illustrate the spatiotemporal evolution of a scroll wave when the L-type Ca 2z conductance has been reduced; to start with, we have a single rotating scroll wave with a straight filament; as time passes, for low Dh (say 10 0 ), the filament of the scroll wave does not break up but it drifts away from the center of the simulation domain; finally it hits the boundaries where it is absorbed. For intermediate to high values of Dh (say 30 0 {60 0 ), the filament of the scroll wave meanders, breaks and leads to a state characterized by scroll-wave turbulence; finally, however, the broken scroll waves are absorbed at the boundaries and the simulation domain is free of scroll waves.  evolution of V (x,y,z,t) is shown, for the representative case with Dh~30 0 , in the video files, Video S4 and Video S5 of the supplementary material; a comparison of Video S4 with Video S3, in which G CaL has its standard value, shows clearly the difference between scroll waves that break up and those that do not. The behaviors of the scroll-wave-filaments is illustrated in figure 5. For  Figure 7 illustrates the spatiotemporal evolution of the tip of the spiral wave, which is obtained by taking a section at z~1 mm through our simulation domain, for the case with Dh~0 0 . In the first sub-figure, namely 7(a.1), which uses the original value of G CaL , the meandering of this tip is limited to a region with linear dimension 0:7 cm, as is expected in a weak-meander régime; in the second of these sub-figures, namely 7(a.2), which uses the modified value of G CaL , the tip meanders considerably as is to be expected in the régime of strong meander. Hence we refer to the respective parameter régimes as the weak meander régime and the strong meander régime in the remaining part of our paper. In most of our simulations, we use a simulation domain of size 9cm|9cm|2mm. We have carried out a few representative studies on a simulation domain of size 13:5cm|13:5cm|2mm that contains 600|600|20 grid points. We find that, in cases where broken scroll waves are eventually absorbed by boundaries, these waves survive longer in large simulation domains than they do in small ones.
If we include transmural heterogeneity in the medium, i.e., our simulation region has epicardial, mid-myocardial, and endocardial domains as we have described in the section entitled Methods, the initially straight filament of the scroll wave is gradually distorted: it bends and twists and finally breaks off at the boundary between the epicardium and the mid-myocardium. We show this in the representative filament-tracking plots in figure 8.

Inhomogeneities
We turn now to a systematic study of scroll-wave dynamics in the 3D TNNP model with FR and inhomogeneities; we consider both conduction and ionic inhomogeneities. As we show below, scroll waves can now exhibit rich dynamical behaviors.

Conduction inhomogeneities
We consider cylindrical conduction inhomogeneities, which we obtain by setting all components of the diffusion tensor equal to a small number (^10 {6 ) in a cylindrical region whose axis is parallel to the z direction and that spans the thickness of our simulation domain. If the radius r of this cylindrical obstacle is small (r 0:5625 cm) the scroll wave separates into two parts when it impinges on the obstacle; these parts sweep around it, rejoin each other beyond it, and emerge finally as a plane wave; thus, the scroll wave does not get anchored to such thin, cylindrical obstacles. Fiber rotation works against such anchoring because it makes different scroll-wave layers, perpendicular to the z axis, rotate at different phases and move with different conduction velocities (CV s): if in one layer the spiral-wave projection of the scroll wave gets anchored successfully to the inhomogeneity, then the spiral-wave projections in other layers are pulled away from it. In two-dimensional planar sections of the scroll wave we see spiral waves; in some layers of the simulation domain these rotate together, whereas in others they rotate out of phase with the wave in some given layer. Furthermore, if we view the scroll wave in three dimensions, we find that its filament distorts in such a way that it remains attached to the obstacle at more than one point along its length in the transmural direction; the spiral sections of such a scroll wave appear pinned to the obstacle in some layers and detached from it in other layers.
Qualitatively different dynamical behaviors are obtained if we increase the radius r of the cylindrical conduction inhomogeneity. To illustrate such behaviors we now consider inhomogeneities with r~1:125 cm. We find that the dynamics of the scroll wave depends sensitively on the position and size of the obstacle, as in our earlier studies without FR [2,3]; furthermore, this dynamics also depends sensitively on the FR angle Dh: If we fix the position of the obstacle, then the dynamical response of the scroll wave depends on Dh; if we fix Dh, then the dynamical response of the scroll wave depends on the position of the obstacle.
Earlier studies of spiral-wave dynamics in two-dimensional mathematical models for cardiac tissue [36] have distinguished between weak-meander and strong-meander régimes. In the former, the filament of the scroll wave (or the tip of the spiral wave in a 2D simulation [3]) does not wander significantly in the simulation domain; in the latter, this filament or tip wandering is significant.
In the sections that follow, we present our results for weak-and strong-meander régimes separately. The parameters used in the original TNNP model [3,44] lead to weak meandering; when we reduce G CaL to one-fourth of its maximal value we obtain the strong-meander régime in the TNNP model.

Weak-meander régime
When the obstacle is located at the center of the domain, the scroll wave anchors to it for Dh~10 0 ,30 0 , and 60 0 , as shown in the representative volume renderings of V (x,y,z,t) in figures 9(a.1), (b.1), and (c.1), respectively; the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S6 of the supplementary material; the corresponding videos for the cases with Dh~10 0 ,30 0 , and 60 0 are available at http:// www.physics.iisc.ernet.in/*rahul/new_movies.html.
If the obstacle lies at a corner, far from the core of the scroll, we find that the scroll wave does not get anchored to the inhomogeneity, as shown in the representative volume renderings of V (x,y,z,t) in figures 9(a.2), (b.2), and (c.2) for Dh~10 0 ,30 0 , and 60 0 , respectively; the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S7 of the supplementary material; the corresponding videos for the cases with Dh~10 0 ,30 0 , and 60 0 are available athttp://www.physics. iisc.ernet.in/*rahul/new_movies.html. If we shift the obstacle position to a corner that is close to the core of the scroll wave, distinct wave-pinning occurs at Dh~10 0 and 30 0 , as illustrated in the representative volume renderings of V (x,y,z,t) in figures 9(a.3) and 9(b.3); the spatiotemporal evolution of V(x,y,z,t) for the representative case of Dh~30 0 is shown in Video S8 of the supplementary material. However, if Dh~60 0 , the delicate balance between anchoring at the inhomogeneity and the FRinduced disruption of wave-pinning is upset: some sections of the wave readily anchor to the inhomogeneity whereas others sweep past the obstacle, reassemble beyond it, and emerge as plane waves; this is shown in the representative volume rendering of V (x,y,z,t) in figure 9(c.3). The spatiotemporal evolution of V (x,y,z,t) for the cases with Dh~10 0 ,30 0 , and 60 0 are shown in   . Representative tip-trajectories of the spiral wave that is obtained by taking a section at z~1 mm through our simulation domain, for the case with Dh~0. In (a.1), which uses G CaL~0 :000175 nS/pF, the meandering of the tip is limited to a region with linear dimension^0:7 cm, as is expected in a weak-meander régime; in (a.2), which uses G CaL~0 :000044 nS/pF, the tip meanders considerably as is to be expected in the régime of strong meander. doi:10.1371/journal.pone.0018052.g007 videos, available at http://www.physics.iisc.ernet.in/*rahul/ new_movies.html.
Spatiotemporal evolution of the filament of the scroll wave in the presence of a cylindrical, conduction-type inhomogeneity located at a corner of the simulation domain, which is initially far from the core of the scroll wave, is illustrated in figure 10. For  Representative plots, at time t~400 and 800 ms of the filament of a scroll wave in the 3D TNNP with transmural heterogeneity and fiber rotation; at t~0 we start with a scroll wave whose filament is straight; the widths of the epicardial, mid-myocardial, and endocardial layers are, respectively, 2:025 mm, 7:2 mm, and 2:025 mm. doi:10.1371/journal.pone.0018052.g008 Figure 9. Scroll-wave dynamics in the 3D TNNP model with fiber rotation and a cylindrical conduction inhomogeneity with radius r~1:125 cm. Obstacle at the center of the domain: the scroll wave anchors to it for Dh~10 0 ,30 0 , and 60 0 , as shown in the representative volume renderings of V (x,y,z,t) in (a.1), (b.1), and (c.1), respectively (the spatiotemporal evolution of V(x,y,z,t) for the representative case of Dh~30 0 is shown in Video S6 in the supplementary material; the videos for the cases with Dh~10 0 ,30 0 and 60 0 are available at http://www.physics.iisc.ernet.in/ *rahul/new_movies.html). Obstacle at a corner, far from the core of the scroll: the scroll wave does not get anchored to the inhomogeneity, as shown in the representative volume renderings of V (x,y,z,t) in (a.2), (b.2), and (c.2) for Dh~10 0 ,30 0 , and 60 0 , respectively (the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S7 in the supplementary material); the videos for the cases with Dh~10 0 ,30 0 and 60 0 are available at http://www.physics.iisc.ernet.in/*rahul/new_movies.html). Obstacle at a corner that is close to the core of the scroll wave: distinct wave-pinning occurs at Dh~10 0 and 30 0 , as illustrated in the representative volume renderings of V(x,y,z,t) in (a.3) and (b.3). However, if Dh~60 0 , some sections of the wave readily anchor to the inhomogeneity whereas others sweep past it, reassemble beyond it, and emerge as plane waves as shown in the volume rendering of V (x,y,z,t) in (c.3) (the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S8 in the supplementary material); the videos for the cases with Dh~10 0 ,30 0 and 60 0 are available at http://www. physics.iisc.ernet.in/*rahul/new_movies.html). doi:10.1371/journal.pone.0018052.g009 dynamics of the scroll wave. The filament, that is located away from the obstacle in each case, displays bending and twisting, but no breakage. Figure 11 illustrates the spatiotemporal evolution of the filament of the scroll wave in the presence of a cylindrical, conduction-type obstacle, located at the center of the simulation domain; for Dh~10 0 , figures 11 (a.1), (a.2), (a.3), (a.4) and (a. 5) show, respectively, the filament of the scroll wave at t~2:4s, t~2:8s, t~3:2s, t~3:6s and t~4s. The analogs of these figures for Initially it remains roughly straight as shown in 11(a.1), but with time, it exhibits bending in the transmural direction, without twisting or breaking up. When Dh~30 0 , the filament anchors to the obstacle. With time, it bends murally as well as transmurally, twisting around the obstacle but without breaking up. When Dh~60 0 , filament breakage occurs in the transmural direction. The broken fragments of the filament remain pinned to the obstacle, but twist around it over and over again. Scroll-wave break-up does not occur away from the obstacle, which explains why our time series data do not show chaotic signatures. Figure 12 illustrates the dynamics of the filament of the scroll wave when the cylindrical, conduction-type obstacle is located at a corner of the simulation domain that was initially close to its core. For Dh~10 0 , figures 12(a.1), and (a.2) show, respectively, the scroll-wavefilament at t~2:4s, and t~3:6s; the analogs of these figures for Dh~30 0 are figures 12 (b.1) and (b.2), and for Dh~60 0 are figures 12(c.1) and (c.2). When Dh~10 0 , the filament is observed to anchor to the obstacle initially; with time it detaches from the obstacle over a major part of its length. A very small fragment is seen to break off from an end of the filament. When Dh~30 0 , the filament is anchored to the obstacle in some layers but not in others. With time, it breaks up transmurally into smaller filaments, which in turn remain attached to the obstacle. When Dh~60 0 , the filament remains detached from the obstacle over a major portion of its length. Gradually with time, breakage occurs in the transmural direction; the broken fragments of the filament get pinned to the obstacle. Scroll-wave break-up does not occur away from the obstacle, which explains why our time-series data do not show chaotic signatures.
If we now include transmural heterogeneity, as described above, and also place a cylindrical conduction inhomogeneity at the center of the simulation domain, a scroll wave, with a filament that is straight initially, breaks up completely in the transmural direction. The wave still remains anchored to the obstacle along its length at various points; but the differences in wave-propagation speeds in the epicardial, mid-myocardial, and endocardial layers lead to distortion and transmural-breakage of the scroll-wave filament. A representative filament plot for this case is shown in figure 13(a.1).

Strong-meander régime
We now turn to studies of scroll-wave dynamics in the strongmeander régime of the 3D TNNP model with conduction inhomogeneities. The value of G CaL is reduced to one-fourth of its maximal value; and the cylindrical conduction-type inhomogeneity has a radius r~1:125 cm.
We begin with the obstacle at the center of the domain. If Dh~10 0 , scroll-wave break up occurs close to the obstacle, broken elements of the scroll wave get pinned to the obstacle, whereas parts of the scroll wave that lie away from the obstacle get annihilated by colliding with each other or they move out of the simulation domain. A small fragment of the scroll wave, which we refer to as a seed wave, breaks off from a region near the core and eventually regenerates a scroll wave; this cycle of the near elimination and subsequent regeneration of the scroll wave is then repeated over and over again. This is illustrated in the representative volume rendering of V (x,y,z,t) in figure 14(a.1) at the instant of time at which the seed wave forms. Figure 14(b.1) displays a volume-rendering plot of the simulation domain (for Dh~30 0 ), at the instant of seed-wave formation during the spatiotemporal evolution of the scroll wave in the presence of a cylindrical, conduction-type heterogeneity at the center of the domain. We find that, if Dh~30 0 , some layers of the domain show sections of the scroll wave anchoring to the obstacle; in other layers the section of the scroll wave is eliminated. Thus, filamentbreakage occurs prominently in the transmural direction. The spatiotemporal evolution of V (x,y,z,t) for this case is shown in Video S9 of the supplementary material and also at http://www. physics.iisc.ernet.in/*rahul/new_movies.html. If we increase Dh to 60 0 , pinning still occurs but for a shorter duration of time than in the cases Dh~30 0 and 10 0 . Parts of the scroll wave are eliminated from most of the layers; this leads to the break up of the scroll-wave filament along the transmural direction, in addition to the break up in the plane; in a few layers a seed wave regenerates a single rotating scroll wave and this process is repeated again and again. We illustrate this in the representative volume rendering of V (x,y,z,t) in figure 14(c.1) at the time instant at which the seed wave forms; the spatiotemporal evolution of V (x,y,z,t) for the representative cases with Dh~10 0 ,30 0 and 60 0 are available at http://www.physics.iisc.ernet.in/*rahul/new_movies.html.
We now place the inhomogeneity near a corner, but far away from the core of the scroll wave. The scroll wave does not get anchored to the inhomogeneity; instead, it impinges on the obstacle, swirls around it, and then is absorbed by the boundaries of the simulation domain. We show this in the representative volume renderings of V (x,y,z,t) in figures 14(a.2), (b.2), and (c.2) for Dh~10 0 ,30 0 , and 60 0 , respectively; the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S10 in the supplementary material; the videos for the cases with Dh~10 0 ,30 0 and 60 0 are available at http://www. physics.iisc.ernet.in/*rahul/new_movies.html.
If we locate the obstacle at a corner that is close to the core of the scroll wave, then we find that for Dh~10 0 the wave leaves the simulation domain completely. If Dh~30 0 , initially the scroll-wave activity is turbulent but the final state consists of a single rotating scroll wave in a fraction of the medium; the rest of the medium remains free of scroll-wave activity ( figure 14b.3)). If Dh~60 0 , the initial state is a turbulent one with broken-scroll-wave elements; the final state comprises a single rotating scroll wave that wanes periodically and is then regenerated by a seed wave. We show these dynamical behaviors in the representative volume renderings

Ionic inhomogeneities
We model cylindrical ionic-type inhomogeneities by setting G CaL~0 :000017 in a cylindrical region of radius r; the axis of this cylinder is parallel to the z direction; and it spans the thickness of our simulation domain. We present results for r 1:125 cm.
In the presence of such an ionic-type inhomogeneity we find that the scroll wave anchors to the obstacle for Dh~10 0 ,30 0 , and 60 0 , when the obstacle is positioned at the center of the simulation domain. We show these dynamical behaviors in the representative volume renderings of V (x,y,z,t) in figures 16 (a.1), (b.1),and (c.1) for Dh~10 0 ,30 0 , and 60 0 , respectively; the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S12 in the supplementary material; the videos for . When Dh~10 0 , the filament is anchored to the obstacle initially; as time passes it detaches from the obstacle over a considerable fraction of its length. A very small fragment is seen to break off from an end of the filament. When Dh~30 0 , the filament is anchored to the obstacle in some layers but not in others. With time, it breaks up transmurally into smaller filaments, which in turn remain attached to the obstacle. When Dh~60 0 , the filament remains detached from the obstacle over a considerable fraction of its length. With the passage of time, breakage occurs in the transmural direction; the broken fragments of the filament get pinned to the obstacle. Scroll-wave break-up does not occur away from the obstacle, which explains why our time-series data do not show chaotic signatures. doi:10.1371/journal.pone.0018052.g012 the cases with Dh~10 0 ,30 0 and 60 0 are available at http://www. physics.iisc.ernet.in/*rahul/new_movies.html. These figures and videos show that, in the presence of the ionic inhomogeneity and FR, the dynamical evolution of the scroll wave is such that the system displays spatiotemporal chaos [2,3] near the location of the inhomogeneity; the rest of the domain supports a scroll wave Figure 13. Scroll-wave filament in the 3D TNNP model with transmural heterogeneity and fiber rotation and cylindrical inhomogeneities. Representative plots, at time t~200ms of the filament of a scroll wave in the 3D TNNP with transmural heterogeneity and fiber rotation; at t~0 we start with a scroll wave whose filament is straight; the widths of the epicardial, mid-myocardial, and endocardial layers are, respectively, 2:025 mm, 7:2 mm, and 2:025 mm; the cylindrical inhomogeneities have radii and heights of 1:125 cm and are located at the center of the simulation domain; left and right panels show conduction and ionic inhomogeneities, respectively. doi:10.1371/journal.pone.0018052.g013 Figure 14. Scroll-wave dynamics in the 3D TNNP model with fiber rotation and a cylindrical conduction inhomogeneity with radius r~1:125 cm and the L-type Ca 2z conductance G CaL reduced to 0:000044 nS/pF from its normal value 0:000175 nS/pF. Obstacle at the center of the domain: (a.1) A representative volume rendering of V (x,y,z,t), with Dh~10 0 , showing broken elements of the scroll wave pinned to the obstacle and the formation of a seed wave. (b.1) A representative volume rendering of V (x,y,z,t), with Dh~30 0 , showing some sections of the scroll wave anchoring to the obstacle in some layers but leaving the domain in other layers and a seed wave (the spatiotemporal evolution of V(x,y,z,t) for this case is shown in Video S9 of the supplementary material, which is also available at http://www.physics.iisc.ernet.in/*rahul/new_movies.html). (c.1) A representative volume rendering of V (x,y,z,t), with Dh~60 0 , showing wave-pinning and eventual wave-elimination from most of the layers and a seed wave regenerating a single rotating scroll wave. Obstacle at a corner, far from the core of the scroll: the scroll wave impinges on the inhomogeneity, swirls around it, and is then absorbed by the boundaries of the simulation domain as shown in (a.2), (b.2), and (c.2) for Dh~10 0 ,30 0 , and 60 0 , respectively (the spatiotemporal evolution of V(x,y,z,t) for the representative case of Dh~30 0 is shown in Video S10 in the supplementary material). Obstacle at a corner that is close to the core of the scroll wave: the wave leaves the simulation domain completely when Dh~10 0 as shown in (a.3); if Dh~30 0 , the final state consists of a single rotating scroll wave in a fraction of the medium and the rest of the medium remains free of scroll-wave activity as illustrated in (b.3); if Dh~60 0 , the final state comprises a single rotating scroll wave that wanes periodically and is then regenerated by a seed wave as shown in (c.3) (the spatiotemporal evolution of V(x,y,z,t) for the representative case of Dh~30 0 is shown in Video S11 in the supplementary material; the videos for each of the cases illustrated in this figure, are available at http://www.physics.iisc.ernet.in/*rahul/ new_movies.html). doi:10.1371/journal.pone.0018052.g014 whose temporal evolution is quasiperiodic. Figure 17 illustrates the spatiotemporal evolution of the filament of the scroll wave, in the presence of a cylindrical, ionic-type obstacle, located at a corner of the simulation domain that was initially far from the core of the scroll wave; for  figure 17(a.3), at a later time, a second filament is formed within the obstacle; the medium then supports more than one scroll wave. When Dh~30 0 , scroll-wave activity is initially supported over a finite region of the domain. The rest of the domain contains plane waves. Within the region that supports scroll-wave activity, the filament is anchored to the obstacle in some layers but not in others. With time, this filament grows in length, detaches from the obstacle and drifts away from it. At an even later time, the filament breaks up transmurally, enters the obstacle and remains pinned to it. When Dh~60 0 , the filament remains detached from the obstacle at all times. Though the filament is initially intact, with time filament breakage occurs in the transmural direction; the broken fragments of the filament also remain detached from the obstacle.
The filament of the scroll in the presence of a cylindrical, ionictype obstacle, located at the center of the simulation domain, is illustrated in figure 18; for When Dh~10 0 , the filament is close to the obstacle initially; with time it attaches to the obstacle and remains pinned to it throughout the duration of the simulation. When Dh~30 0 , the filament initially twists around the obstacle, without attaching to it.
With time, this filament enters the obstacle and successfully anchors to it. When Dh~60 0 , the filament remains twisted around the obstacle at all times, without actually anchoring on to it. Breakage occurs close to the obstacle, in the transmural direction. Spatiotemporal evolution of the filament of the scroll wave in the presence of a cylindrical, ionic-type inhomogeneity, located at a corner of the simulation domain, which was initially close to the core of the scroll wave, is illustrated in figure 19; for Dh~10 0 , figures 19(a.1), (a.2), and (a.3) show, respectively, the scroll-wave filament at t~2:4s, t~3:6s, and t~4s. The analogs of these figures for Dh~30 0 , are figures 19(b.1), (b.2), and (b.3), and for Dh~60 0 are figures 19(c.1), (c.2), and (c.3). When Dh~10 0 , the filament initially anchors to the obstacle, briefly detaches from it and then reanchors to the obstacle. In the process, a second scroll hook forms within the inhomogeneity, which exists side by side with the original filament. When Dh~30 0 , the filament of the scroll wave is anchored to the obstacle in some layers but not in others. The formation of a figure-eight-like structure occurs within the inhomogeneity; this later disappears and the medium acquires the potential to support more than one scroll; hence the existence of three filaments is observed. When Dh~60 0 , the filament initially bends in the transmural direction, without twisting. Then it enters the obstacle, breaks up and forms a scroll ring which eventually moves out of the inhomogeneity and merges with the main filament to form a long filament that twists angularly around the obstacle.
Space  (x~45 mm,y~45 mm,z~1mm), which lies inside the inhomogeneity; the disorder in the space-time plots and the broad-band background in the power spectra, both signatures of spatiotemporal chaos, increase with Dh. If we record data from a point outside the ionic inhomogeneity then we find no indication of spatiotemporal chaos as can be seen from the illustrative plots, for Dh~30 0 , in figures 21(a.1) and 21(a.2); the former is a spacetime (z,t) pseudocolor plot of V (x,y,z,t) from the representative point (x~22:5mm,y~67:5mm) and the latter a power spectrum E(v) obtained from the time series V (x,y,z,t) from the representative point (x~22:5mm,y~67:5mm,z~1mm); the peaks in this power spectrum can be indexed as n 1 v 1f zn 2 v 2f , with v 1f^3 :5 and v 2f^7 :5 two incommensurate frequencies, so the temporal evolution of the scroll wave, away from the ionic inhomogeneity, is quasi-periodic. This coexistence of spatiotemporal chaos and quasi-periodic behavior in different parts of the simulation domain has been observed earlier in numerical studies of the 2D TNNP model with ionic inhomogeneities [3] but never before in studies of 3D mathematical models for cardiac tissue.
If the ionic inhomogeneity is located near a corner that is away from the core of the scroll wave, for low Dh 30 0 , the scroll wave does not get anchored to the inhomogeneity; instead, it impinges on the obstacle, goes around it, reassembles as a plane wave beyond it, and, far from the inhomogeneity, it evolves as it did in the absence of the inhomogeneity. Figure 21(b.1) shows a spacetime (z,t) pseudocolor plot of V (x,y,z,t) from the representative point (x~22:5mm,y~67:5mm) and 21(b.2) is a power spectrum E(v) obtained from the time series V (x,y,z,t) from the representative point (x~22:5mm,y~67:5mm,z~1mm); when the ionic inhomogeneity is located near a corner, which is away from the core of the scroll wave, and Dh~30 0 . However, if Dh~60 0 , different layers of the simulation domain support sections of the scroll wave that rotate quasi-periodically with different conduction velocities. We show this in the representative for Dh~10 0 ,30 0 , and 60 0 , respectively; the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S13 in the supplementary material; the videos for the cases with Dh~10 0 ,30 0 and 60 0 are available at http://www. physics.iisc.ernet.in/* rahul/new_movies.html. If we place the cylindrical ionic inhomogeneity near a corner of the simulation, which is also close to the core of the scroll wave, we find that the wave gets pinned to the inhomogeneity for all the values of Dh we consider, namely, 10 0 ,30 0 , and 60 0 ; in the case Dh~60 0 the scroll wave leaves the simulation domain in some layers but continues rotating in other layers. We show this in the representative volume renderings of V (x,y,z,t) in figures 16(a.3), (b.3), and (c.3) for Dh~10 0 ,30 0 , and 60 0 , respectively; the spatiotemporal evolution of V (x,y,z,t) for the representative case of Dh~30 0 is shown in Video S10 in the supplementary material; the videos for the cases with Dh~10 0 ,30 0 and 60 0 are available at http://www.physics.iisc.ernet. in/*rahul/new_movies.html. Figure 21(c.1) shows the space-time (z,t) pseudocolor plot of V (x,y,z,t) from the representative point (x~22:5mm,y~67:5mm) and 21(c.2) shows the power spectrum E(v) obtained from the time-series V (x,y,z,t) from the representative point (x~22:5mm,y~67:5mm,z~1mm); when the ionic inhomogeneity is located near a corner, which is also close to the core of the scroll wave, and Dh~30 0 .
If we now include transmural heterogeneity, as described above, and also place a cylindrical ionic inhomogeneity at the center of the simulation domain, a scroll wave, with a filament that is straight initially, breaks up transmurally at the boundary between the epicardium and the mid-myocardium. The scroll-wave filament enters the region occupied by the ionic inhomogeneity, leaves it, moves away from it, but always remains close to it. A representative filament plot for this case is shown in figure 19(a.2).

Discussion
We have presented a systematic study of the combined effects of muscle-fiber rotation and inhomogeneities on scroll-wave dynamics in the TNNP model [44] for cardiac tissue. This model, which has been introduced recently, is based on experimental data obtained from human ventricular cells; it allows for variations of intracellular ion concentrations, contains 12 ionic currents, 12 gating variables, one ion pump, and an ion exchanger; all major ionic currents are included here, e.g., the fast inward Na z current, the L-type Ca 2z current, the transient-outward potassium current, the slow-potassium-delayed-rectifier current, the rapid-potassiumdelayed-rectifier current, and the inward-rectifier K z current as described in detail in Ref. [3]. To the best of our knowledge, no other study has attempted to systematize scroll-wave dynamics, in the presence of fiber rotation and inhomogeneities, in such a realistic mathematical model for cardiac tissue. The quantitative results from our study have been described in detail in the previous section. Here we discuss the important qualitative issues that emerge from our study.
We have shown in earlier work [2,3] that spiral-wave dynamics depends sensitively on the position, shape, and size of inhomogeneities in mathematical models for cardiac tissue; these studies have been restricted to two-dimensional simulation domains, for realistic models like the TNNP model; in three-dimensional simulation domains only illustrative studies have been carried out with inhomogeneities for the simple, two-variable Panfilov model [3,54]. Our study here extends considerably the results of Refs. [2,3] by using the three-dimensional TNNP model with fiber rotation. We find, in this case, that, in addition to a sensitive dependence on the positions, shapes, sizes, and types of inhomogeneities, scroll-wave dynamics also depends delicately upon the degree of fiber rotation, which is measured in our calculations by the angle Dh.
In one set of studies, we have used parameters, such as the values of certain channel conductances, that are appropriate for the epicardium [44]; whereas, in another set, which contains three representative simulations, we have used parameters that are appropriate for the epicardium, the mid-myocardium as well as the endocardium [44].
In the first set of studies, we have used a slab thickness of 2 mm because the approximate width of the epicardium in the left ventricle of a human heart is 2 mm (see, e.g., [51]). We have then carried out a study that has never been attempted before, namely, we have elucidated the combined effects of fiber rotation and inhomogeneities on scroll-wave dynamics in such tissue. Studies such as those of Refs. [20,32,36] have concentrated on the effects of fiber rotation on scroll-wave dynamics but without inhomogeneities; our earlier work has elucidated the effects of different types of inhomogeneities -conduction and ionic -on spiral-and scrollwave dynamics in mathematical models of cardiac tissue but without fiber rotation. To the best of our knowledge, the study we present here is the first one to investigate the combined effects of fiber rotation and inhomogeneities on scroll-wave dynamics in the detailed ionic TNNP model (principally with epicardial parameters). In the second set of studies, we have included transmural heterogeneity. Hence we have used a simulation domain of total thickness 1:125 cm.
We have checked in earlier studies (see Ref. [3]) that most qualitative aspects of spiral-wave dynamics in the presence of a conduction inhomogeneity do not depend sensitively on whether we use Neumann boundary conditions on the inhomogeneity or whether we merely set D equal to a very small number inside this inhomogeneity; this method of reducing D has also been used in Ref. [52]. In this study we are not investigating in detail the interaction of the tip of the spiral wave with the inhomogeneity and the removal of the spiral from such an inhomogeneity by pacing; this might well require a smooth gradient in D as suggested in Refs. [50,53]. . When Dh~10 0 , the filament is free from the obstacle initially; it remains unaffected by the presence of the obstacle throughout the duration of the simulation. However, as shown in figure (a.3), at a later time, a second filament is formed within the obstacle; the medium then supports more than one scroll wave. When Dh~30 0 , scroll-wave activity is initially supported over a finite region of the domain. Within the region that supports scroll-wave activity, the filament anchors to the obstacle in some layers but not in others. With the passage of time, this filament grows in length, detaches from the obstacle and drifts away from it. At an even later time, the filament breaks up transmurally, enters the obstacle, and remains pinned to it. When Dh~60 0 , the filament remains detached from the obstacle at all times; initially it is straight; with the passage of time breakage occurs in the transmural direction; the broken fragments of the filament also remain detached from the obstacle. doi:10.1371/journal.pone.0018052.g017 We find that scroll waves do not anchor easily to cylindrical conduction inhomogeneities with small radii; however, the tendency for anchoring increases with the radius. If parameters such as G CaL have values that lead to a weak meandering of the scroll wave, then there is a fine balance between anchoring of this wave at the inhomogeneity and a disruption of wave-pinning by fiber rotation; this disruption increases with Dh. If the parameters are such that the system is in the strong-meander régime, then again the anchoring is suppressed as we increase Dh; furthermore, the scroll wave can be eliminated from most of the layers only to be regenerated by a seed wave over and over again. Ionic inhomogeneities can also lead to an anchoring of the scroll wave; again fiber rotation works against this anchoring tendency. Scroll waves cannot enter the region inside a conduction inhomogeneity but they can enter the region inside an ionic inhomogeneity; thus, in the latter case, we obtain far richer scroll-wave dynamics than in the former. In particular, we can get a coexistence of spatiotemporal chaos and quasi-periodic behavior in different parts of the simulation domain; this has not been noted before in studies of 3D mathematical models for cardiac tissue. In the simple, conventional picture persistent arrhythmias arise because of the presence of heterogeneities that act as pinning centers; thus, cardiac tissue with a high density of large, in-excitable obstacles is more likely to exhibit VT than tissue that does not. This simple picture may have to be revised in the light of our study.
The formation of a seed wave takes place only when three conditions are simultaneously satisfied: (a) the model parameters are such that the system lies in the strong-meander regime; (b) fiber rotation is present; and (c) the scroll wave interacts with an obstacle whose radius is large enough to anchor it. Given these conditions, if we look near the filament of the scroll wave, we find that the waveback is not smooth: it develops bumps whose 2D analogs have been noted for the LRI model in Refs. [16,36]. In this 2D study the conduction velocity near the tip varies between fast and slow phases; our 3D study exhibits a similar phenomenon near the scroll-wave filament as it meanders. When it encounters an obstacle, the filament impinges on it and tries to anchor there but the oscillations of the conduction velocity near the core are significantly large; the slow parts of the wave (representing the excitations in the part of the tissue that is relatively refractory) are snapped off the parent wave, while it is attached to the obstacle; this part that has been snapped off is the seed wave.
Major advances in biomedical engineering have found substantial applications in the clinical diagnosis and treatment of cardiac arrhythmias and the electro-physiological study of waves of electrical activation in cardiac tissue and their role in arrhythmias. However, the methods that are used to record such waves in experiments on cardiac tissue are typically constrained to probe only the surface of the tissue. Our in silico study of scroll waves in the realistic 3D TNNP model has allowed us to elucidate, by detailed three-dimensional visualization, the effects of inhomogeneities and fiber rotation on the dynamics of these waves. Such visualization is not easy in experiments on three-dimensional cardiac-tissue preparations. Thus, our studies offer insights that complement those derived from experimental measurements [17,[55][56][57][58][59][60] of spiral and scroll waves in cardiac tissue.
The development of a detailed understanding of the influence of fiber rotation and inhomogeneities on scroll-wave activity in cardiac tissue has at least two clinical implications. Such an understanding is important in developing efficient, low-amplitude defibrillation schemes [2,3,61,62]. Furthermore, clinicians should eventually be able to assess the risk of cardiac arrhythmias in different patients from fiber-orientation data, obtained from diffusion-tensor magnetic-resonance imaging (DTMRI).
There are a few limitations of our study: We do not use a bi-domain model, and we do not use an anatomically realistic simulation domain. We conjecture that neither one of these limitations affect the qualitative results of our study; in particular, we expect that, even if we use a bidomain version of the TNNP model and an anatomically realistic simulation domain, the sensitive dependence of scroll-wave dynamics on inhomogeneities and fiber rotation will persist. Also, recent work [63] has shown that the differences between mono-domain and bidomain mathematical models for cardiac tissue are not significant in the absence of large electrical stimuli.
Representative in silico studies of transmural heterogeneities in 2D, without the types of obstacles we concentrate on and without fiber rotation, can be found in Refs. [41][42][43]. The inclusion of transmural heterogeneity, in three-dimensional simulation domains with fiber rotation, is of importance in developing a detailed understanding of scroll-wave dynamics in ventricular tissue; it has been studied explicitly and systematically for the TNNP model in Ref. [35] but without the types of inhomogeneities we consider here. To perform a detailed, systematic, study of scroll-wave dynamics in the TNNP model, with fiber rotation, transmural heterogeneity and different types of inhomogeneities, such as the ones we have given above, lies beyond the scope of this paper.
However, we present three representative results from our simulations on the 3D TNNP model with transmural heterogeneity, fiber rotation, both with and without conduction and ionic inhomogeneities; such studies are computationally expensive.

Methods
We use the TNNP model [44] of human cardiac tissue for our in silico studies. This model is defined by the following reactiondiffusion equation for the transmembrane potential V : here I ion is the total ionic current density, which is expressed as a sum of six major and six minor ionic currents: I major~INa zI CaL zI to zI Ks zI Kr zI K1 I minor~INaCa zI NaK zI pCa zI pK zI bNa zI bCa ; . When Dh~10 0 , the filament initially anchors to the obstacle, briefly detaches from it, and then anchors again to the obstacle. In the process, a second scroll hook forms within the inhomogeneity, which exists alongside the original filament. When Dh~30 0 , the filament of the scroll wave anchors to the obstacle in some layers but not in others; within the inhomogeneity the filament forms a figure-eight-type structure, which later disappears; the medium now shows more than one scroll-wave filament. When Dh~60 0 , the filament initially bends in the transmural direction, without twisting. Then it enters the obstacle, breaks up, and forms a scroll ring which eventually moves out of the inhomogeneity and merges with the main filament to form a long filament that twists around the obstacle. doi:10.1371/journal.pone.0018052.g019 here I Na is the fast inward Na z current, I CaL the L{type slowinward Ca 2z current, I to the transient outward current, I Ks the slow, delayed-rectifier current, I Kr the rapid, delayed-rectifier current, I K1 the inward rectifier K z current, I NaCa the Na z =Ca 2z exchanger current, I NaK the Na z =K z pump current, I pCa the plateau Ca 2z current, I pK the plateau K z currents, I bNa the background Na z current, and I bCa the background Ca 2z current. Here time t is in milliseconds, voltage V in millivolts, current densities I X in picoamperes per picofarad (pA/pF), conductances (G X ) in nanoSiemens per picofarad (nS/pF), and the intracellular and extracellular ionic concentrations (X i , X o ) in millimoles per liter (mM/L). When we write I X , we implicitly mean the total current density per unit capacitance density, because we measure all ionic currents as current density per unit capacitance per unit area as in second-generation models (see, e.g., Refs. [44,[64][65][66]). The capacitance density is measured in units of pF/cm 2 , so the unit of I X , as we have used it, is (pA/cm 2 ) per (pF/cm 2 ), which is nothing but pA/pF. Area and capacitance are related because the specific capacitance of cardiac tissue is of the order of 1mF/cm 2 . For a detailed list of the parameters of this model and the equations that govern the spatiotemporal behaviors of the transmembrane potential and currents here, we refer the reader to Refs. [3,44]. In our numerical simulations of the TNNP model we use a threedimensional domain that is a rectangular slab of dimension 9 cm |9 cm |2 mm. We use a finite-difference method in space with grid spacings Dx~Dy~0:0225 cm and Dz~0:01 cm and a seven-point stencil for the Laplacian; for time marching we use an Euler method with a time step Dt~0:02ms. We model cardiac-muscle-fiber orientation as in Refs. [32,67]: The simulation domain is viewed as a stack of planes with fibers arranged parallel to each other in any given plane. The fibers in each plane are rotated, with respect to those in adjacent planes, by an angle a; thus, the orientation of the fibers in a plane is specified by an angle h(z) that varies continuously with the transmural coordinate z. In a normal human heart the average total FR, from the endocardium to the epicardium, is^110 0 {120 0 . The distribution of this angle over the three-layered heart wall varies across mammalian species and, within a specie, from individual to individual. In our study we consider, in most cases, only a patch of tissue from the epicardium, so we use a smaller value of this average total FR, which we denote by Dh. We vary Dh to analyze the effects of FR on the dynamics of scroll waves in the TNNP model for cardiac tissue, both with and without inhomogeneities. We incorporate the effects of FR in this model by using a diffusion tensor D of the form used in Refs. [32,36]; this has the following components: Here D 11~DE cos 2 h(z)zD \1 sin 2 h(z),D 22~DE sin 2 h(z)z D \1 cos 2 h(z), and D 12~D21~(  (b.1) shows the same for the case where the obstacle is at a corner away from the core of the scroll; and (c.1) shows the same for the case where the obstacle is at a corner that is close to the scroll-wave core. There is no indication of spatiotemporal chaos as can be seen from these illustrative space-time plots. Power spectra E(v) obtained from the time series V(x,y,z,t) recorded from the representative point (x~22:5mm,y~67:5mm,z~1mm) for these three cases are illustrated in (a.2), (b.2), and (c.2); the peaks in these power spectra can be indexed as n 1 v 1f zn 2 v 2f , with v 1f^3 :5 and v 2f^7 :5 two incommensurate frequencies, so the temporal evolution of the scroll wave, away from the ionic inhomogeneity, is quasi-periodic. doi:10.1371/journal.pone.0018052.g021 the diffusivity for propagation parallel to the fiber axis, D \1 the diffusivity perpendicular to this axis but in the same plane, and D \2 the diffusivity perpendicular to the fiber axis but out of the plane, i.e., in the transmural direction. h(z), the twist angle along the transmural direction, can expressed as h(z)~al z , where a is the rate of fiber rotation and l z the distance of a given layer from the bottommost layer of the slab of tissue. We assume that h(z) changes continuously in such a way that the total fiber rotation across the tissue is Dh. In our studies we use Dh~0 0 ,10 0 ,30 0 and 60 0 . We use D E~0 :00154 cm 2 /ms and D \1~D\2~0 :0002 cm 2 /ms.
The initial condition is generated by stacking spiral waves, obtained from 2D simulations on the TNNP model [3], one on top of the other along the z axis, so as to generate a simple scroll wave with a straight filament. The system is then allowed to evolve in time (i) in the absence of and (ii) in the presence of localized heterogeneities.We then record the spatiotemporal evolution of the transmembrane potential V .
We consider two types of inhomogeneities: (a) the first is a conduction inhomogeneity; (b) and the second is an ionic inhomogeneity. The former can arise because of an inactive clump of dead cells or scar tissue. Ionic inhomogeneities can arise from functional blockage and the presence of deposits. We model case (a) by introducing a solid, isopotential cylinder of radius r and voltage {86:2 mV at different positions in the simulation domain; the axis of the cylinder is chosen to be along the z direction; and all elements of the diffusion tensor are assigned very small values (^10 {6 cm 2 /ms) within the cylinder. The radii of these cylinders are varied. In our simulations we have used r~0:5625cm to model small obstacles and r~1:125 cm to model large obstacles. To model the ionic inhomogeneities in case (b) we modify the value of G CaL in similar, but not isopotential, cylindrical regions.
We can also include the effects of transmural heterogeneity in the TNNP model. This requires a thicker simulation domain than the one we have used above. The typical thickness of the human ventricular wall is 1{2 cm and the total fiber rotation Dh^110{120 0 . This wall is divided into three layers, namely, the epicardium, the mid-myocardium, and the endocardium. The boundary between the epicardium and the mid-myocardium is sharp, but that between the endocardium and the mid-myocardium is irregular; in our simulations we do not account for such irregularity. We divide our simulation domain into three slabs placed one on top of the other; the thicknesses of the epicardial, mid-myocardial, and endocardial slabs are chosen to be 2:025 mm, 7:2 mm, and 2:025 mm, respectively; for each one of these slabs we use parameters, such as the transient outward channel conductance G to , that are appropriate for epicardial, mid-myocardial, and endocardial regions; the values of these parameters are given in Ref. [44]. We have performed three representative simulations with such transmural heterogeneity: (a) in the first we do not have any localized inhomogeneities; (b) in the second we include a cylindrical, conduction-type inhomogeneity; and (c) in the third we introduce a cylindrical, ionic-type inhomogeneity; for cases (b) and (c) the cylindrical inhomogeneities, with radii and heights of 1:125 cm, are placed at the center of the simulation domain. Our simulations with transmural heterogeneity show more complex intramural filament shapes and broken filaments than those we observe for scroll-wave filaments in thin epicardial simulation domains. Transmural heterogeneity is perhaps the principal cause for this increase in complexity; but the increase in the thickness of the simulation domain and the degree of fiber rotation are also important. Detailed simulations, which lie beyond the scope of this paper, will have to be designed to see whether the effects of transmural heterogeneity, tissue thickness, and fiber rotation on scroll-wave dynamics can be disentangled from each other. Furthermore, in order to obtain scroll-wave break up in the TNNP model, with fiber rotation (and/or inhomogeneities), we reduce the Calcium channel conductance, G CaL to one-fourth of its original value. This, in effect, blocks I CaL . Several experimental and computational studies have shown that altering the value of G CaL can induce transitions from weak to strong meandering of the tip of a spiral wave as, e.g., in the studies of Refs. [16,36] for the LRI model. A transition from a single rotating spiral to spiralwave break up can be obtained by altering other types of ionic conductances also. For example, in our previous studies [3] of the two-dimensional TNNP model of ventricular tissue we have shown that increasing the plateau Calcium conductance G pCa or decreasing the plateau Potassium conductance G pK can be responsible for such spiral-wave break-up. Here we change G CaL and the total fiber rotation (FR) and study their effects on scrollwave dynamics in our three-dimensional simulation domain because the Calcium ion channel can be controlled easily by pharmaceutical means as has been done in several studies [68][69][70].

Supporting Information
Video S1 Scroll-wave dynamics in the 3D TNNP model in the absence of fiber rotation and inhomogeneities: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 1  Video S6 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation Dh~30 0 and a cylindrical conduction inhomogeneity of radius r~1:125 cm, located at the center of the simulation domain: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 5 (b.1), with 2ƒtƒ4s and 15 frames per second.

(MPEG)
Video S7 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 and a cylindrical conduction inhomogeneity of radius r~1:125 cm, located at a corner of the simulation domain, far from the core of the scroll: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 5 (b.2), with 2ƒtƒ4s and 15 frames per second.

(MPEG)
Video S8 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 and a cylindrical conduction inhomogeneity of radius r~1:125 cm located at a corner of the simulation domain that is close to the scroll-wave core: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 5  Video S9 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 and a cylindrical conduction inhomogeneity of radius r~1:125 cm located at the center of the simulation domain, and with the L{type Ca 2z conductance G CaL , reduced to 0:000044 nS/pF from its normal value 0:000175 nS/pF: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 6 (b.1), with 2ƒtƒ4s and 15 frames per second.

(MPEG)
Video S10 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 and a cylindrical conduction inhomogeneity of radius r~1:125 cm located at the corner of the simulation domain, far from the core of the scroll, and with the L{type Ca 2z conductance G CaL , reduced to 0:000044 nS/pF from its normal value of 0:000175 nS/pF: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t) for the same parameter values as in figure 6 (b.2), from 2ƒtƒ4s and with 15 frames per second. (MPEG) Video S11 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 , a cylindrical conduction inhomogeneity of radius r~1:125 cm located at the corner of the simulation domain that is close to the scroll-wave core, and with the L{type Ca 2z conductance G CaL , reduced to 0:000044 nS/ pF from its normal value 0:000175 nS/pF: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 6 (b.3), with 2ƒtƒ4s and 15 frames per second. (MPEG) Video S12 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 and a cylindrical ionic inhomogeneity of radius r~1:125 cm located at the center of the simulation domain: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 8 (b.1), with 2ƒtƒ4s and 15 frames per second. (MPEG) Video S13 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 and a cylindrical ionic inhomogeneity of radius r~1:125 cm located at a corner of the simulation domain far from the core of the scroll: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 8  Video S14 Scroll-wave dynamics in the 3D TNNP model in the presence of fiber rotation with Dh~30 0 and a cylindrical ionic inhomogeneity of radius r~1:125 cm located at a corner of the simulation domain that is close to the scroll-wave core: An animated volume rendering, illustrating the spatiotemporal evolution of the transmembrane potential V (x,y,z,t), for the same parameter values as in figure 8 (b.3), with 2ƒtƒ4s and 15 frames per second. (MPEG)