advertisement

Nonequilibrium Arrhythmic States and Transitions in a Mathematical Model for Diffuse Fibrosis in Human Cardiac Tissue Rupamanjari Majumder1, Alok Ranjan Nayak1, Rahul Pandit1,2* 1 Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore, India, 2 Jawaharlal Nehru Centre for Advanced Scientific Research, Bangalore, India Abstract We present a comprehensive numerical study of spiral-and scroll-wave dynamics in a state-of-the-art mathematical model for human ventricular tissue with fiber rotation, transmural heterogeneity, myocytes, and fibroblasts. Our mathematical model introduces fibroblasts randomly, to mimic diffuse fibrosis, in the ten Tusscher-Noble-Noble-Panfilov (TNNP) model for human ventricular tissue; the passive fibroblasts in our model do not exhibit an action potential in the absence of coupling with myocytes; and we allow for a coupling between nearby myocytes and fibroblasts. Our study of a single myocytefibroblast (MF) composite, with a single myocyte coupled to Nf fibroblasts via a gap-junctional conductance Ggap , reveals five qualitatively different responses for this composite. Our investigations of two-dimensional domains with a random distribution of fibroblasts in a myocyte background reveal that, as the percentage Pf of fibroblasts increases, the conduction velocity of a plane wave decreases until there is conduction failure. If we consider spiral-wave dynamics in such a medium we find, in two dimensions, a variety of nonequilibrium states, temporally periodic, quasiperiodic, chaotic, and quiescent, and an intricate sequence of transitions between them; we also study the analogous sequence of transitions for three-dimensional scroll waves in a three-dimensional version of our mathematical model that includes both fiber rotation and transmural heterogeneity. We thus elucidate random-fibrosis-induced nonequilibrium transitions, which lead to conduction block for spiral waves in two dimensions and scroll waves in three dimensions. We explore possible experimental implications of our mathematical and numerical studies for plane-, spiral-, and scroll-wave dynamics in cardiac tissue with fibrosis. Citation: Majumder R, Nayak AR, Pandit R (2012) Nonequilibrium Arrhythmic States and Transitions in a Mathematical Model for Diffuse Fibrosis in Human Cardiac Tissue. PLoS ONE 7(10): e45040. doi:10.1371/journal.pone.0045040 Editor: Xiongwen Chen, Temple University, United States of America Received June 18, 2012; Accepted August 11, 2012; Published October 8, 2012 Copyright: ß 2012 Majumder et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: The study was funded by the Council for Scientific and Industrial Research, India, and the Department of Science and Technology, India. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. * E-mail: rahul@physics.iisc.ernet.in inhomogeneities [1,5–21]. Some of these studies [5,6,14–18] concentrate on the interaction of a spiral or scroll wave with a localized inhomogeneity; others are devoted to investigations of the effects of a large number of randomly-distributed, point-type inexcitable obstacles on such waves [1,7–13]. We concentrate on the latter types of studies here because they have been designed to mimic arrays of fibrotic strands or diffuse fibrosis in cardiac tissue. By using a simple model for cardiac tissue with many inexcitable obstacles, Pertsov [7] has shown that such obstacles can influence the rotation of spiral waves and lead to anisotropies in propagation. Turner, et al. [8] have studied the effects of fibrosis in the Priebe-Beukelmann model. Spach, et al. [9] have used the Nygren model for human atrial tissue and mimicked the effects of diffuse fibrosis by removing lateral gap junctions; they find that with such heterogeneity in intercellular couplings, there is a tendency for partial wave block and re-entry. Kuijpers, et al. [10] have used the Courtemanche model for human atrial tissue and heterogeneous uncoupling to model diffuse fibrosis. These studies indicate that fibrosis can increase vulnerability to re-entry; however, they have not explored in detail the effects of fibrosis on the dynamics of spiral and scroll waves in these models. Such Introduction Extra-cellular-matrix (ECM) materials constitute about 6% of the volume of human cardiac tissue in an average, healthy heart [1]. These include fibroblasts, non-excitable collagen, and elastin fibrils, which fill the subepicardial space between the epicardium and myocardium [2] and bridge the gaps between myocardial tissue layers. The major component of the ECM are fibroblast cells that produce interstitial collagen, of types I, III, IV, and VI [3]. These contribute to myocardial structure, cardiac development, cell-signaling, and electro-mechanical functions in myocardial tissue. In mammalian cardiac tissue, fibroblast cells show an intimate spatial interrelation with every myocyte that borders one or more fibroblasts [4]. In tissue containing both myocytes and fibroblasts, it has been assumed traditionally that gap-junctional couplings exist exclusively between myocytes; but recent experimental studies have shown that there is a functional, heterogeneous, myocyte-fibroblast coupling [3,4]. Computer simulations of electrical-wave propagation in mathematical models for cardiac tissue have been used to investigate the interplay of spiral and scroll waves with conduction and other PLOS ONE | www.plosone.org 1 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions oscillatory state in which the initial AP response to the external stimulus is followed by oscillations of the membrane potential about a mean value without the application of any other external stimulus. In régime R5, the MF composite produces a single AP under the influence of the external stimulus; after that it does not return to the resting state but to another time-independent state in which it is non-excitable. We then study propagation of plane waves in a 2D simulation domain with TNNP-type [16,23] myocytes (M) or fibroblasts (F) of the type described in Ref. [24]; M and F are distributed randomly through the simulation domain; and there are diffusive couplings between nearest-neighbor cells. We investigate plane-wave propagation through both mural slices, with epicardial parameters, and transmural slices, consisting of epicardial, mid-myocardial, and endocardial regions, in moderate- and strong-coupling cases, i.e., with myocyte-to-fibroblast diffusion constants of 0:0000218cm2 =ms and 0:000048cm2 =ms, respectively. We obtain stability diagrams for both these cases in the Pf {Ef parameter space. In the moderate-coupling case, this stability diagram is simple: for low values of Pf the plane wave leaves the system, which returns to an excitable state; for large values of Pf the plane wave is annihilated by target waves and the medium is left in a state that is weakly excitable or not excitable at all. In the strong-coupling case the stability diagram is very rich; it contains the spatiotemporal analogs of the régimes R1–R5 mentioned above for an isolated MF composite. The last part of our study examines the effect of diffuse fibrosis on spiral-wave dynamics in 2D and scroll-wave dynamics in 3D with myocytes and fibroblasts distributed randomly as above; we concentrate on the strong-coupling case here. At low values of Pf , we find that single, rotating spiral and scroll waves have slightly corrugated wave fronts, but they propagate much as they do in the absence of fibroblasts. For large values of Pf , we find that such spiral and scroll waves do not propagate through the simulation domain and are either (a) annihilated by spontaneously generated target waves or (b) absorbed at the boundaries. This crossover from a state with propagating waves and electrical activity to a state with no electrical activity occurs via a sequence of nonequilibrium transitions; the precise sequence depends on the initial conditions and the realization of the disordered array of M and F cells. Given the spatial and temporal resolution we have been able to achieve in 2D, we find the following rough sequence of states: at low Pf we begin with a state with a single spiral rotating periodically (SRSP); as Pf increases this gives way to a state with a single spiral rotating quasi-periodically (SRSQ); as Pf increases we obtain a state with multiple spirals that rotate periodically (MRSP); this then gives way to a state with a multiple spirals that rotate quasi-periodically (MRSQ), which is followed by a spiral-turbulence (ST) state and, eventually, by the absorption state (SA). In 3D, the analogous sequence of states, which we have been able to resolve, is the following: at low Pf we begin with a state with a single rotating scroll wave (SRS); as Pf increases this gives way to a state with multiple rotating scroll waves (MRS); this then gives way to the absorption state (SA). an exploration has been initiated by Panfilov [11] and ten Tusscher, et al. [1,12,13], who investigate the effects of diffuse fibrosis on the propagation of electrical waves of activation and arrhythmogenesis in both two-variable and detailed ionic mathematical models for human ventricular tissue; they model fibrosis as non-conducting inhomogeneities that are distributed randomly in their simulation domain. They show that, as the concentration of such inhomogeneities increases, CV decreases for plane-wave propagation, the wave fronts become jagged, and there is an increase in the tendency for the formation and break up of spiral waves; at sufficiently large densities of these inhomogeneities, they find that complete conduction blockage occurs. McDowell, et al. [22] have used a three-dimensional computational model, based on MRI data, of chronically infarcted rabbit ventricles to characterize arrhythmogenesis because of myofibroblast infiltration as a function of myofibroblast density; this study includes periinfarct zones (PZ), ionic-current remodeling therein, and different degrees of myofibroblast infiltration. Their work shows that, at low densities, myofibroblasts do not alter the propensity for arrhythmia; at intermediate densities, myofibroblasts cause AP shortening and thus increase this propensity; at high densities, these myofibroblasts protect against arrhythmia by causing resting depolarization and blocking propagation. We present a major extension of the work of ten Tusscher, et al. [1,12,13] on diffuse fibrosis in mathematical models for cardiac tissue by introducing fibroblasts randomly in the state-of-the-art TNNP model [16,23] for human ventricular tissue; the fibroblasts are passive, insofar as they do not exhibit an action potential in the absence of coupling with myocytes. Our model for the fibroblasts is much more realistic than the one used by ten Tusscher, et al. [1,12,13]; in particular, we use the fibroblast model of MacCannell, et al. [24,25]; we also allow coupling between nearby myocytes and fibroblasts. The parameters in this model cannot be determined precisely from experiments [4,25] so it is important to explore a wide, but biophysically relevant, range of parameters. Our in silico study is well suited for such an exploration so it is very effective in complementing experimental studies of electrical-wave propagation in fibrotic cardiac tissue. We begin with an overview of our principal results before we present the details of our study. We first study a single myocytefibroblast (MF) composite in which a single myocyte is coupled to Nf fibroblasts via a gap-junctional conductance Ggap . We study two cases, namely, moderate and strong coupling between fibroblasts and myocytes; for each one of these cases, we consider three parameter sets [23] for the myocytes that are suitable for epicardial, mid-myocardial, and endocardial layers of the heart wall; experiments suggest that 0:1nSƒGgap ƒ8nS [26,27]; thus, we consider Ggap ~0:1 nS, Ggap ~4 nS, and Ggap ~8 nS to be the weak-, moderate- and strong-coupling cases, respectively. We excite each such MF composite via an electrical stimulus and then record its responses for different values of Nf and Ef . In the moderate-coupling case (Ggap ~4 nS), the electrical load of the fibroblasts on the myocyte is not significant, except in a very narrow range of parameters. However, in the strong-coupling case (Ggap ~8 nS), for different ranges of the parameters Nf and Ef , we observe five qualitatively different responses for the MF composite; we call them R1–R5. In R1 the MF composite responds essentially like an uncoupled myocyte. In régime R2, the MF composite produces a secondary AP, after the first one that is generated by the external stimulus. In R3, this composite displays autorhythmicity, i.e., it fires a train of APs, after the first external stimulus, and without the application of subsequent stimuli; each AP in this autorhythmic train differs from the normal AP of an uncoupled mycocyte. In régime R4, the MF composite displays an PLOS ONE | www.plosone.org Materials and Methods The first system we study is a single myocyte-fibroblast (MF) composite in which a single myocyte is coupled to Nf fibroblasts via a gap-junctional conductance Ggap ; we consider the range 1ƒNf ƒ10. We then carry out studies in 2D and 3D simulation domains in which myocytes M or fibroblasts F are distributed randomly through the simulation domain; we include diffusive 2 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions millivolts, conductances (GX ) in nanoSiemens per picofarad (nS/ pF), the intracellular and extracellular ionic concentrations (Xi , Xo ) in millimoles per liter (mM/L) and current densities, per unit capacitance, IX in picoamperes per picofarad (pA/pF), as used in second-generation models (see, e.g., Refs. [23,29–31]). For a detailed list of the parameters of this model and the equations that govern the spatiotemporal behaviors of the transmembrane potential and currents, we refer the reader to Refs. [16,23]. For the fibroblasts, we use the model of MacCannell, et al. [24]; i.e., we treat the fibroblasts as passive circuit elements that couple with the myocyte in the MF composite; and the fibroblast ionic i is current Iion,f couplings between these. The precise ionic models and the diffusive couplings are described in detail below. In 2D we use a square simulation domain (13:5cm|13:5cm), when we consider a mural slice, and a rectangular domain (13:5cm|1:35cm), when we study a transmural slice. This 1:35cm thick transmural slice is further divided into three parallel strips of width 2:7mm for the epicardium, 7:875mm for the midmyocardium, and 2:925mm for the endocardium. For the mural slice, we choose parameters for epicardial-type myocytes, whereas, for the transmural slice, we use parameters for epicardial, midmyocardial, and endocardial myocytes in the appropriate strips. In 3D our simulation domain is a rectangular parallelepiped of physical size 9:36cm|9:36cm|1:125cm. Our spatial grid uses dx~dy~0:0225 cm in 2D and dx~dy~dz~0:0225 cm in 3D; and our time-marching scheme has a time step dt~0:02 ms. The actual thickness of the human left ventricular wall varies between 1 and 2 cm (see Ref. [28]); we have chosen a thickness of 1.35 cm because it is well within this range. For our 3D simulations, we have reduced this thickness to 1.125 cm for the following reasons: (a) we have checked, in a few representative cases, that the principal results of our study do not change qualitatively if we reduce the thickness of the domain by a few millimeters; (b) our 3D simulations are computationally very expensive because we need to run the simulations for long time and store several intermediate configurations to distinguish the states SRS, MRS, and SA; (c) a wall thickness of 1.125 cm also lies within the acceptable range of thicknesses for the left ventricular wall. The precise ratios of the thicknesses of the three layers of the heart wall are not known. However, a rough estimate (see Ref. [28]) indicates that the epicardium is, on average, 2–3 mm thick, the mid-myocardium is the thickest zone, and the endocardium has a highly non-uniform thickness; the values we have chosen for the thicknesses of the epicardial, mid-myocardial, and endocardial slices in our simulations are commensurate with these rough estimates and observations (see Ref. [17]). For the ionic activity of the myocytes, we adopt the state-of-theart TNNP model [23] for human cardiac tissue, which is the following reaction-diffusion equation for the transmembrane potential V : LV Iion ~+:(D+V ){ ; Lt Cm i ~Gf (Vfi {Ef ), Iion,f where Gf , Vfi and Ef , are, respectively, the conductance, transmembrane potential, and the resting membrane potential for the fibroblast. We incorporate muscle-fiber anisotropy in both 2D and 3D simulations, as in Refs. [32,33]; we account for diffusive couplings between nearest-neighbor myocytes, nearest-neighbor fibroblasts, next-nearest neighbor myocytes, next-nearest-neighbor fibroblasts and nearest-neighbor myocyte-fibroblast pairs. We use two diffusion tensors, namely, Dmm and Dff , for myocyte-myocyte (mm) and fibroblast-fibroblast (ff ) diffusive couplings. The diffusion tensors Dmm and Dff have the form used in Refs. [32,34]; we give this form below for a diffusion tensor that is denoted generically by D and has, in three dimensions, the components shown hereunder: 2 D11 6 D~4 D21 0 3 0 7 0 5, D\2 D11 ~DE cos2 h(z)zD\1 sin2 h(z), ð1Þ D22 ~DE sin2 h(z)zD\1 cos2 h(z), D12 ~D21 ~(DE {D\1 ) cos h(z) sin h(z): ð4Þ For myocyte-myocyte couplings, we use D~Dmm ; mm 2 DE ~0:00154 cm =ms is the diffusivity for propagation through 2 myocytes, parallel to the fiber axis, Dmm \1 ~0:000385 cm =ms the diffusivity through myocytes, perpendicular to this axis but in the 2 same plane, and Dmm \2 ~0:000385 cm =ms the diffusivity through myocytes, perpendicular to the fiber axis but out of the plane, i.e., in the transmural direction. The twist angle along the transmural direction is h(z). It is related to the rate of fiber rotation a via h(z)~alz , where lz is the thickness of the tissue as measured from the bottom of the endocardium, along the z axis. The total fiber rotation (FR) across the slab is taken to be 1100 , which is the typical value for the human ventricular wall. For fibroblastfibroblast couplings Dff also has the same form as D; we choose DffE ~0:000385 cm2 =ms and Dff\1 ~Dff\2 ~0:00009265 cm2 =ms because we expect ff diffusive couplings to be much smaller than their mm counterparts. In 2D simulations, we use only the ð2Þ Imajor ~INa zICaL zIto zIKs zIKr zIK1 ; Iminor ~INaCa zINaK zIpCa zIpK zIbNa zIbCa ; here INa is the fast inward Naz current, ICaL the L{type slowinward Ca2z current, Ito the transient outward current, IKs the slow, delayed-rectifier current, IKr the rapid, delayed-rectifier current, IK1 the inward rectifier K z current, INaCa the Naz =Ca2z exchanger current, INaK the Naz =K z pump current, IpCa the plateau Ca2z current, IpK the plateau K z currents, IbNa the background Naz current, and IbCa the background Ca2z current. The time t is measured in milliseconds, voltage V in PLOS ONE | www.plosone.org D12 D22 0 where here Iion , the total ionic current, is expressed as a sum of the following six major and six minor ionic currents: Iion ~Imajor zIminor ; ð3Þ 3 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions components D11 ,D12 ,D21 , and D22 for a mural slice in the x{y plane; for a transmural slice in the x{z plane we retain only the components D11 and D\2 . We turn now to mf and fm diffusive couplings; the magnitudes of these are not known well experimentally, nor has the role of fiber orientation been investigated in this context. Therefore, in the interests of a parsimonious description, we neglect fiber orientation in the mf and fm diffusive couplings Dmf and Dfm , respectively, and treat them as scalars. In keeping with the idea that the interaction between a myocyte and a fibroblast should be weaker than that between two myocytes, but stronger than that between two fibroblasts, we use the following illustrative values: (a) for the strong-coupling case Dfm ~0:00141cm2 =ms and Dmf ~0:000048cm2 =ms; and (b) for the moderate-coupling case Dfm ~0:000642cm2 =ms and Dmf ~0:0000218cm2 =ms; in both these cases Dfm ~Dmf (Cm =Cf ), where the total cellular capacitances for myocytes and fibroblasts are Cm ~185 pF and Cf ~6:3 pF, respectively [24]. Note that, in our 2D and 3D models, there is no on-site coupling between myocytes and fibroblasts; this has been translated into diffusive couplings between such cells if they are at nearest-neighbor sites in our simulation domains. We generate the initial condition for our studies by using the following protocols: We begin with only myocytes on all sites of our 2D simulation domain. For plane-wave-propagation studies, we apply a stimulus, of amplitude 150pA=pF for 2 ms, along one edge of the simulation domain. For our studies of spiral-wave dynamics, we obtain a spiral wave in the 2D domain by using the method proposed by Shajahan, et al. [16]. In our 3D scroll-wave studies, we begin with an initial scroll wave that consists of our initial, 2D spiral waves stacked one on top of the other; thus, we begin with a simple scroll wave with a straight filament as in Ref. [17]. Pseudocolor plots of Vm for these spiral and scroll waves, which we use as initial conditions in our subsequent studies, are given in Figs. 1 (a) and (b), respectively. For every value of Pf , we generate a random array of myocytes and fibroblasts in our 2D and 3D simulation domains as follows by using a random-number generator to assign F or M to a site such that the percentage of F sites is Pf ; illustrative arrays of F and M are given in Fig. 2 (a)–(e). This distribution of F and M cells is held fixed throughout the subsequent spatiotemporal evolution of the initial spiral and scroll waves described above; in the language of condensed-matter physics, a static configuration of F and M is an example of quenched disorder [35–37]. At the initial time, the fibroblast transmembrane potential Vf is set equal to its resting value Ef ~{49:6 mV [24]. The temporal evolution of the transmembrane potential VA of the cell at site A in the lattice is governed by LVA {Iion,A ~ zDA , Lt CA ð5Þ where DA indicates the diffusion term. This can be written most easily in discrete form and it depends on (a) whether the cell at site A is M or F and (b) whether the cell on the neighboring site is of type M or F. We illustrate the form of the diffusion term D for a two-dimensional mural slice for three representative sites A, B, and C in Fig. 2(f), which have M, M, and F cells, respectively: DA ~D11 mm ( zDmf ( DC ~Dfm ( (dx)2 )zDmm 22 ( V (i,jz1){V (i,j) (dy)2 ð6Þ ð6Þ ) V (i,j){V (i,j{1) (V (iz1,jz1){V (iz1,j{1)) ; )zD12 mm ( 2dxdy (dy)2 DB ~Dmf ( zDmf ( V (iz1,j)zV (i{1,j){2V (i,j) V (iz1,j)zV (i{1,j){2V (i,j) (dx)2 )zDmm 22 ( V (i,jz1){V (i,j) (dy)2 ð7Þ ð7Þ ) V (i,j){V (i,j{1) (V (i{1,j{1){V (i{1,jz1)) ; )zD12 mm ( 2dxdy (dy)2 V (iz1,j)zV (i{1,j){2V (i,j) (dx)2 )zDff22 ( V (i,jz1)zV (i,j{1){2V (i,j) (dy)2 ): ð8Þ ð8Þ In our studies with the MF composite, we apply a stimulus current of 52pA=pF for 1 ms to the composite and allow the system to evolve in time. We record the membrane potential of the central myocyte in the MF composite and plot it as a function of time for different values of Nf and Ef . For our measurements of CV and l, we prepare the 2D simulation domain as discussed above and initiate a plane wave by applying a stimulus of amplitude 150pA=pF for 2 ms along the left edge (y axis) of the domain. We record the time series of Vm at four representative points of the domain. For studies on the mural slice, these points are (3:375cm,3:375cm), (10:125cm,3:375cm), (3:375cm,10:125cm) and (6:75cm,6:75cm); for studies on the transmural slice, these points are (3:375cm,0:3375cm), (3:375cm,1:0125cm), (10:125cm,0:3375cm), (10:125cm,1:0125cm) and (6:75cm,0:675cm). From these time-series data, we obtain the times t1 and t2 at which the upstroke of the action potential (AP) is initiated at pairs of points that are separated by Dx along the axis parallel to the direction of propagation of the wave; CV ~Dx=Dt, where Dt~t2 {t1 ; the wavelength l~CV APD90% , where APD90% is the action-potential duration at 90% repolarization; we obtain average values for CV and l over the four points mentioned above. Results In this Section we present the results of our computational studies. We begin with our investigation of MF composites and discuss how their action potential is influenced by the number of fibroblasts Nf , their resting membrane potential Ef , the gapjunctional coupling Ggap , and the myocyte parameters, which distinguish myocytes from the epicardium, the mid-myocardium, and the endocardium. We then explore the propagation of plane waves of electrical activation through 2D simulation domains with randomly distributed myocytes and fibroblasts such that the percentage of fibroblasts is Pf ; we consider propagation through Figure 1. The initial configurations for the spiral and scroll waves in our 2D and 3D simulations (see text). doi:10.1371/journal.pone.0045040.g001 PLOS ONE | www.plosone.org 4 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 2. Spatial distributions of myocytes and fibroblasts in our simulation domain. Illustrative plots of the distributions of fibroblasts and myocytes, on a representative 60|60 part of our simulation domain in two dimensions (2D). Fibroblasts (filled black squares) and myocytes (unfilled white squares) are distributed randomly in our simulation domain so that a given grid point contains either a myocyte or a fibroblast, but not both; the percentage of fibroblasts Pf is (A) 5%, (B) 15%, (C) 25%, (D) 35%, and (E) 45%. In (F) we show an illustrative arrangement of myocytes (pink circles) and fibroblasts (blue ellipses); for such an arrangement the diffusion terms are as given in Eqs. 6–8. doi:10.1371/journal.pone.0045040.g002 both mural and transmural slices. Next we consider spiral-wave dynamics in 2D and 3D simulation domains with Pf % fibroblasts. The action potential durations (APDs) are different, for uncoupled myocytes from the endocardium, the mid-myocardium and the epicardium. The APD of the myocyte-fibroblast composite (MF) is also different for the three types of myocytes. However, the APD of an MF depends not only on the type of myocyte but also on the values of the gap-junctional conductance Ggap , the resting membrane potential of fibroblasts (Ef ), and the number of fibroblasts coupled to a myocyte (Nf ). Fig. 3 shows pseudocolor plots of the APD for MFs, with epicardial, mid-myocardial, and endocardial myocytes, as functions of Ggap and Ef , for Nf = 1, 2, 3, and 4. In our studies, Ggap is moderate (4 nS) or strong (8 nS), so the influence of the gap-junctional coupling is quite significant. When Nf = 1, the differences between the APDs, for the three types of MFs, is considerable, for all values of Ggap and Ef ; as Nf increases, this difference between the APDs is significant only at low values of Ggap ; in particular, for Nf = 4 and Ggap = 8 nS, the distinction between these APDs is almost negligible. However, if Nf = 4 and Ggap v2, this distinction is quite prominent at all the values of Ef that we have considered. For studies on transmural heterogeneity in mouse tissue, please refer to [38,39]. PLOS ONE | www.plosone.org MF Composite The response of the myocyte-fibroblast (MF) composite depends on Nf ,Ef , and Ggap and the properties of the myocyte. In both moderate- and strong-coupling cases, with Ggap ~4 nS and Ggap ~8 nS, respectively, if Nf ~1 the MF composite produces a single action potential on the application of an external stimulus and then returns to the normal resting membrane potential for myocytes (^{86mV); we designate this as a response of type R1; and we illustrate this, for the case Ggap ~8 nS, by plots of Vm versus time t in Figs. 4 (a.1), (a.2), and (a.3) for epicardial, midmyocardial, and endocardial myocytes, respectively. Four other types of responses are possible and are listed below and portrayed in Fig. 4, for the case Ggap ~8 nS: R2: In this case there is a secondary AP, after the first one generated by an external stimulus; the MF composite then returns to the resting state as in R1 (Figs. 4 (b.1), (b.2), and (b.3) for epicardial, mid-myocardial, and endocardial myocytes, respectively). R3: Here the MF composite is autorhythmic, i.e., it produces a sequence of APs, after the first external stimulus; each AP in this autorhythmic sequence differs from the normal AP of an uncoupled mycocyte (Figs. 4 (c.1), (c.2), and (c.3) for epicardial, mid-myocardial, and endocardial myocytes, respectively). R4: The MF composite can display an oscillatory response in which the initial, stimulus-induced AP is 5 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 3. Dependence of action potential duration (APD) of an MF composite on Ggap and Ef . Illustrative pseudocolor plots of the APD for MFs, with epicardial, mid-myocardial, and endocardial myocytes, as functions of the gap-junctional conductance Ggap and the resting membrane potential of fibroblasts Ef , for the number of coupled fibroblasts Nf = 1, 2, 3, and 4. doi:10.1371/journal.pone.0045040.g003 followed by oscillations of Vm , about a mean value, without the application of any other external stimulus (Figs. 4 (d.1), (d.2), and (d.3) for epicardial, mid-myocardial, and endocardial myocytes, respectively). R5: The MF composite produces a single AP because of an external stimulus; after that it does not return to the normal resting state but to another time-independent state in which it is non-excitable (Figs. 4 (e.1), (e.2), and (e.3) for epicardial, mid-myocardial, and endocardial myocytes, respectively). The regions in which our MF composite displays responses of types R1–R5 are shown, for illustrative parameter values, in the Ef {Nf plane in Figs. 5 (a.1), (a.2), (b.1), (b.2), (c.1) and (c.2). For the moderate-coupling case, Ggap ~4 nS, Figs. 5 (a.1), (b.1), and (c.1) show the stability diagrams for, respectively, endocardial, PLOS ONE | www.plosone.org mid-myocardial, and epicardial myocytes in the MF composite; their analogs for the strong-coupling case, Ggap ~8 nS, are given in Figs. 5 (a.2), (b.2), and (c.2); here régimes R1, R2, R3, R4, and R5 are denoted, respectively, by yellow hexagrams, red squares, green circles, pink diamonds, and blue pentagrams. All these régimes appear in stability diagrams for the moderate- and strongcoupling cases; but régimes R3 and R4 occupy very small areas especially in the moderate-coupling case; and régime R5, which occupies a significant fraction of the stability diagram in the strong-coupling cases, occurs in a narrow parameter range in the case of moderate coupling, but only when we consider MF composites with mid-myocardial myocytes. 6 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions PLOS ONE | www.plosone.org 7 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 4. Representative time series of the transmembrane potential recorded from an MF composite. Plots of the transmembrane potential Vm of the myocyte in the myocyte-fibroblast (MF) composite versus time t for the strong-coupling case Ggap ~8 nS showing the following: responses of type R1 for Nf ~1, Ef ~{49 mV, and (a.1) epicardial myocytes, (a.2) mid-myocardial myocytes, and (a.3) endocardial myocytes; responses of type R2 for (b.1) epicardial myocytes and Nf ~4, and Ef ~{21 mV, (b.2) mid-myocardial myocytes and Nf ~5, and Ef ~{32 mV, and (b.3) endocardial myocytes and Nf ~4, and Ef ~{22 mV; autorhythmic responses of type R3 for (c.1) epicardial myocytes and Nf ~5, and Ef ~{34 mV, (c.2) mid-myocardial myocytes and Nf ~2, and Ef ~{20 mV, and (c.3) endocardial myocytes and Nf ~4, and Ef ~{29 mV; oscillatory responses of type R4 for (d.1) epicardial myocytes and Nf ~3, and Ef ~{2 mV, (d.2) mid-myocardial myocytes and Nf ~4, and Ef ~{19 mV, and (d.3) endocardial myocytes and Nf ~3, and Ef ~{9 mV; responses of type R5 for (e.1) epicardial myocytes and Nf ~6, and Ef ~{30 mV, (e.2) mid-myocardial myocytes and Nf ~5, and Ef ~{13 mV, and (e.3) endocardial myocytes and Nf ~8, and Ef ~{16 mV. doi:10.1371/journal.pone.0045040.g004 Figure 5. Stability diagrams in the Ef {Nf parameter space for the responses of an MF composite. The regions in which the MF composite shows responses of types R1–R5 are shown in (a.1), (a.2), (b.1), (b.2), (c.1), and (c.2). For the moderate-coupling case, Ggap ~4 nS, (a.1), (b.1), and (c.1) show the stability diagrams for, respectively, endocardial, mid-myocardial, and epicardial myocytes in the MF composite; their analogs for the strong-coupling case, Ggap ~8 nS, are given in (a.2), (b.2), and (c.2); the régimes R1, R2, R3, R4, and R5 are denoted, respectively, by yellow hexagrams, red squares, green circles, pink diamonds, and blue pentagrams. doi:10.1371/journal.pone.0045040.g005 PLOS ONE | www.plosone.org 8 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Plane-wave propagation in a 2D domain We now investigate the propagation of plane waves of electrical activation through a 2D simulation domain of the type we have described above. In this domain, we distribute myocytes and fibroblasts randomly such that the percentage of fibroblasts is Pf ; we consider propagation through both mural and transmural slices. In addition to Pf , the other important parameters in this part of our study are Ef and the components of the diffusion ff mf tensors, i.e., Dmm and Dfm (see Eqs. 4). Recall that, ij , Dij , and D in our 2D and 3D models, there is no on-site coupling between myocytes and fibroblasts; but we have diffusive couplings between such cells if they are at nearest-neighbor sites in our simulation domains; here Pf plays a rôle similar to that of Nf , in our studies of MF composites. We show below that the temporal responses, of types R1–R5, for MF-composites, have spatiotemporal analogs when we consider plane-wave propagation through our 2D simulation domain; we denote these analogs by the same symbols, namely, R1–R5, because the spatiotemporal evolution of the plane waves in these stability régimes can be rationalized, qualitatively, in terms of the responses of MF composites that we have discussed above. We first consider plane-wave propagation through a mural slice. We find the five qualitatively different spatiotemporal behaviors R1–R5. In the régime R1, which occurs both in moderate- and strong-coupling cases, the plane wave propagates smoothly through the simulation domain but with a slightly corrugated wave front. In the régime R2, small clusters of fibroblasts can form around some sites with myocytes; these clusters have the capacity to generate one subsidiary action potential (cf. the response R2 of an MF composite), before returning to a resting potential that is above the resting potential of the myocytes; because of this subsidiary action potential, target waves are generated by the fibroblast clusters and a plane wave, which tries to propagate through the simulation domain, is annihilated by these target waves, so the whole domain returns to a potential that is above the normal resting membrane potential of myocytes, but below their threshold potential; R2 is absent in the moderate-coupling case, in the parameter régimes that we have explored. The parameter régime R3 is characterized by autorhythmicity; the fibroblast clusters about some myocyte now develop the ability to sustain rhythms of their own (cf. the response R3 of the MF composite); here too the initial plane wave is annihilated by the target waves that are generated by the autorhythmic fibroblast clusters; however, unlike in the case R2, the activity of our medium does not stop here; after a considerable length of time, the autorhythmic fibroblast clusters generate sustained beats of their own; beats from fibroblast clusters of different sizes, which are in different parts of the medium, are incoherent; R3 is absent in the moderatecoupling case, in the parameter régimes that we have explored. In the oscillatory regime régime R4 (cf. the response R4 of the MF composite), the fibroblast clusters produce an initial target wave that annihilates the plane wave; this is followed by temporal oscillations, about some mean potential, of the local membrane potential; R4 is absent in the moderate-coupling case, in the parameter régimes that we have explored. In régime R5 the initial plane wave is terminated by collisions with numerous target waves, which are generated by the fibroblast clusters that are distributed randomly throughout the medium; once the plane wave is removed, the medium moves into a quiescent state with a membrane potential that lies above the excitation-threshold potential for an uncoupled myocyte; no further excitation is possible; R5 is absent in the moderate-coupling case, in the parameter régimes that we have explored. The stability diagram, PLOS ONE | www.plosone.org Figure 6. Stability diagram in the Ef {Pf plane for plane-wave propagation through a mural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. The stability diagram shows the regions with spatiotemporal behaviors R1–R5 in the strong-coupling case (Dmf ~0:000048cm2 =ms); the regions R1, R2, R3, R4, and R5 are denoted, respectively, by blue diamonds, green triangles, pink pentagrams, black squares, and red circles; the spatiotemporal evolution of plane waves in these regions is described in the text. doi:10.1371/journal.pone.0045040.g006 which shows the regions with spatiotemporal behaviors R1–R5 in the strong-coupling case, is given in Fig. 6; regions R1, R2, R3, R4, and R5 are denoted, respectively, by blue diamonds, green triangles, pink pentagrams, black squares, and red circles. Representative pseudocolor plots of the local membrane potential V (x,y,t) are given in Fig. 7 for several values of the time t to illustrate plane-wave propagation, through a 2D mural slice, in the moderate-coupling case, for different values of Pf . (V (x,y,t)~Vm (x,y,t), if the site (x,y) is occupied by a myocyte, and V (x,y,t)~Vf (x,y,t), if the site (x,y) is occupied by a fibroblast.) We do not see behaviors of types R2–R5 here; the plane wave propagates through the medium with a slightly corrugated wave front (region R1). Video S1 shows the spatiotemporal evolution of the plane waves in Figs. 7 (a.1), (b.1), (c.1), (d.1), (e.1), and (f.1). Analogous plots, for the strong-coupling case, of plane-wave propagation, through a 2D mural slice, are shown in Fig. 8; planewave propagation for régime R1 is illustrated in Figs. 8 (a.1)–(a.4), for régime R2 in Figs. 8(b.1)–(b.4), for régime R3 in Figs. 8(c.1)– (c.4), for régime R4 in Figs. 8 (d.1)–(d,4), and for régime R5 in Figs. 8(e.1)–(e.5). The spatiotemporal evolution of these plane waves is given in Video S2. We turn now to illustrative studies of plane-wave propagation through a 2D transmural slice. Here too, we find the five qualitatively different spatiotemporal behaviors R1–R5. In the régime R1, which occurs both in moderate- and strong-coupling cases, the plane wave propagates smoothly through the simulation domain but with remarkable distortion. In the moderate-coupling case, at low values of Pf , the wavefront acquires a smoother appearance than in the 2D mural slice; the smoothness begins to disappear as Pf increases. Furthermore, these waves propagate differently within the three layers of the heart wall, inside the simulation domain. For sufficiently large values of Pf , electrical conduction is partially blocked in the mid-myocardium and completely blocked in the endocardium; the excitation then travels 9 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 7. Pseudocolor plots of the local membrane potential V illustrating plane-wave propagation through a mural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. Here we consider the moderate-coupling case Dmf ~0:0000218 cm2/ms, and Ef ~{30mV and in (a.1)–(a.5) the percentage of fibroblasts Pf ~5%, in (b.1)–(b.5) Pf ~10%, in (c.1)–(c.5) Pf ~15%, in (d.1)–(d.5) Pf ~20%, in (e.1)–(e.5) Pf ~25%, and in (f.1)–(f.5) Pf ~30%. (For full spatiotemporal evolutions see Video S1.). doi:10.1371/journal.pone.0045040.g007 (region R1). Video S3 shows the spatiotemporal evolution of the plane waves in Figs. 9 (a.1), (b.1), (c.1), (d.1), (e.1), and (f.1). Analogous plots, for the strong-coupling case, of plane-wave propagation, through a 2D transmural slice, are shown in Fig. 11; plane-wave propagation for régime R1 is illustrated in Figs. 11 (a.1)–(a.4), for régime R2 in Figs. 11(b.1)–(b.4), for régime R3 in Figs. 11(c.1)–(c.4), for régime R4 in Figs. 11 (d.1)–(d,4), and for régime R5 in Figs. 11(e.1)–(e.5). The spatiotemporal evolution of these plane waves is given in Video S4. only along the epicardium as illustrated in Fig. 9. In the strongcoupling case, régime R1 occurs only at low values of Pf , as in the moderate-coupling case; and here the wave has a smooth wave front. Régimes R2, R3, R4 and R5, analogous to those in the strong-coupling case of the 2D mural slice, are also observed in the strong-coupling case of the 2D transmural slice. However, these are absent in the moderate-coupling case, in the parameter régimes that we have explored. The stability diagram, which shows the regions with spatiotemporal behaviors R1–R5 in the strong-coupling case, is given in Fig. 10; regions R1, R2, R3, R4, and R5 are denoted, respectively, by blue diamonds, green triangles, pink pentagrams, black squares, and red circles. Representative pseudocolor plots of the local membrane potential V (x,y,t) are given in Fig. 9 for several values of the time t to illustrate plane-wave propagation, through a 2D transmural slice, in the moderate-coupling case, for different values of Pf . We do not see behaviors of types R2–R5 here; the plane wave propagates, through the medium, with a distorted wave front PLOS ONE | www.plosone.org Dependence of the conduction velocity and the wavelength on the percentage of fibrosis We characterize the influence of fibroblasts on plane-wave propagation through our mathematical model for myocardial tissue with fibroblasts by studying the dependence of the planewave-conduction velocity CV and the wavelength l on the percentage of fibrosis Pf ; we present illustrative studies at a fixed value of the resting membrane potential of fibroblasts, namely, 10 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 9. Pseudocolor plots of the local membrane potential V illustrating plane-wave propagation through a transmural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. Here we consider the moderate-coupling case Dmf ~0:0000218 cm2/ms, and Ef ~{30mV and in (a.1)–(a.5) the percentage of fibroblasts Pf ~0%, in (b.1)–(b.5) Pf ~5%, in (c.1)–(c.5) Pf ~15%, in (d.1)–(d.5) Pf ~25%, in (e.1)–(e.5) Pf ~35%, and in (f.1)– (f.5) Pf ~40%. As Pf increases, not only does the distortion of the wavefront increase but the wave also propagates preferentially through the zone that has epicardial myocytes (rather than the zones with midmyocardial and endocardial myocytes).(For full spatiotemporal evolutions see Video S3.). doi:10.1371/journal.pone.0045040.g009 Figure 8. Pseudocolor plots of the local membrane potential V illustrating plane-wave propagation through a mural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. Here we consider the strong-coupling case and different régimes in the stability diagram of Fig. 6; (a.1)–(a.4) propagation in régime R1 (Pf ~5%,Ef ~{25mV ); (b.1)–(b.4) propagation in régime R2 (Pf ~10%,Ef ~{25mV); (c.1)–(c.4) propagation in régime R3 (Pf ~15%,Ef ~{25mV); (d.1)–(d.4) propagation in régime R4 (Pf ~20%,Ef ~{25mV ); (e.1)–(e.4) propagation in régime R5 (Pf ~25%,Ef ~{25mV).(For full spatiotemporal evolutions see Video S2.). doi:10.1371/journal.pone.0045040.g008 V (x,y,t), recorded from a point near the corner of the simulation domain, i.e., from (x~2:25cm,y~2:25cm), is shown alongside in Fig. 13(b); Fig. 13(c) shows the power spectrum E(v) of this time series; and the corresponding plot of the inter-beat interval (IBI) versus the beat number n is depicted in Fig. 13(d). The simple periodicity of this time series, the appearance of a single, major peak in E(v) at the fundamental frequency vf ^4Hz, and the constancy of the IBI confirm that the spiral wave in SRSP evolves completely periodically in time. Next we increase Pf in steps of 1%. For Pf *; 14%, the system continues in the state SRSP; but, as Pf approaches 14%, the single, completely periodic, spiral-wave develops a granular texture that increases with Pf ; the distance from the wave-front to the wave-back also decreases. In Fig. 14 we show, for representative values of Pf , pseudocolor plots of the local transmembrane potential V (x,y,t); these plots illustrate the time evolution of a spiral wave in six qualitatively different states, namely, SRSP, SRSQ, MRSP, MRSQ, ST, and SA, which we have defined above; the spatiotemporal evolution of V (x,y,t) for these states is shown in Video S5. The states SRSP and SRSQ have single spirals that rotate periodically and quasiperiodically, respectively; MRSP and MRSQ have multiple spirals whose temporal evolution is periodic and quasiperiodic, respectively; the state ST displays spiral-wave turbulence; and in SA the spiral wave is absorbed at the boundaries of our simulation domain. To examine the temporal evolution of spiral waves in these states, it is useful to look at time series of V (x,y,t), from representative points in the simulation domain, and the resulting plots of the IBI and the power spectra E(v). These are shown for illustrative values of Pf in Figs. 15 and 16. In Figs. 15 (a.1)–(d.3) we have chosen the values of Pf so that we can show examples of temporal 2-cycles (Figs. 15 (a.1)–(a.3) for Pf ~21%), 3-cycles (Figs. 15 (b.1)–(b.3) for Pf ~21:3%), 4-cycles (Figs. 15 (c.1)–(c.3) for Pf ~16:8%), and 5-cycles (Figs. 15 (d.1)– Ef ~{40mV; we choose this value because, from our single-MFcomposite studies, it is evident that, at such a moderately low value of Ef , the MF composite responds to electrical stimuli as in the régime R1, so it is convenient to measure CV . When Pf ~0% we find that CV ^70cm=s, the typical value for plane-wave propagation through the human myocardium; and l^19cm. As we increase Pf , in the moderate-coupling case, CV decreases gradually, as does l. When the MF diffusive coupling is strong, CV decreases gradually at first, but then, once the fibroblast clusters become large enough to generate target waves that can annihilate the plane wave, CV falls rapidly to zero. The medium then may or may not show conduction blockage, depending on whether it has passed into the régime R5, or is still in R2, R3, or R4. Plots of CV and l versus Pf are given, respectively, in Figs. 12 (a) and (b), for both moderate-coupling (open blue circles) and strong-coupling (filled black circles) cases. Influence of diffuse fibrosis on spiral waves in 2D e now explore the dynamics of spiral waves of electrical activation in our mathematical model in the presence of fiber anisotropy and diffuse fibrosis. We start with a monolayer of myocytes (Pf ~0%) and the initial condition of Fig. 1 (a); we observe that, even after t~20 s, the medium supports only one, temporally periodic, rotating spiral wave, which shows no breaks. We call this state SRSP (Single-Rotating-Spiral-Periodic); a representative pseudocolor plot of the local membrane potential V (x,y,t), is given in Fig. 13(a) for t~10 s; the time series of PLOS ONE | www.plosone.org 11 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 10. Stability diagram in the Ef {Pf plane for plane-wave propagation through a transmural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. The stability diagram shows the regions with spatiotemporal behaviors R1– R5 in the strong-coupling case (Dmf ~0:000048cm2 =ms); the regions R1, R2, R3, R4, and R5 are denoted, respectively, by blue diamonds, green triangles, pink pentagrams, black squares, and red circles; the spatiotemporal evolution of plane waves in these regions is described in the text. doi:10.1371/journal.pone.0045040.g010 (d.3) for Pf ~20:3%); these cycles show up most clearly in the IBI plots (Figs. 15 (a.2), (b.2), (c.2), and (d.2)) but their presence can also be surmised from the time series of V (Figs. 15 (a.1), (b.1), (c.1), and (d.1)) and the sharp peaks in the power spectra (Figs. 15 (a.3), (b.3), (c.3), and (d.3)). In Figs. 16 (a.1)–(d.3) we have chosen the values of Pf so that we can show examples of temporal 6-cycles (Figs. 16 (a.1)–(a.3) for Pf ~21:1%), 7-cycles (Figs. 16 (b.1)–(b.3) for Pf ~20:7%), 9-cycles (Figs. 16 (c.1)–(c.3) for Pf ~16:6%), and 10-cycles (Figs. 16 (d.1)– (d.3) for Pf ~16:1%); these cycles show up most clearly in the IBI plots (Figs. 16 (a.2), (b.2), (c.2), and (d.2)) but their presence can also be surmised from the time series of V (Figs. 16 (a.1), (b.1), (c.1), and (d.1)) and the sharp, fundamental frequencies in the power spectra (Figs. 16 (a.3), (b.3), (c.3), and (d.3)). PLOS ONE | www.plosone.org Long time series are required to ascertain the temporal periodicity of these states. Here we obtain local time series for V (x,y,t), from the representative point (x~6:75cm,y~6:75cm), for 0ƒtƒ20 s, which corresponds to 106 time steps; to remove the effects of initial transients, it is best to disregard data from the first 300000 iterations or so. Given plots such as those of Fig. 15 and 16, we can systematize the sequence of transitions that leads from the state SRSP to SA. For the initial conditions and the distributions of fibroblasts that we use, the sequence of transitions is shown in Fig. 17(a). The exact sequence in which these transitions occur depends sensitively on the initial conditions, boundary effects, and the realizations of fibroblast distributions within the domain, as in other nonequilibrium transitions (see, e.g., Refs. [40–42]). 12 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 13. Spiral-wave dynamics in the absence of fibroblasts (Pf ~0%) in a 2D simulation domain with myocytes of the midmyocardial type. (A) Illustrative pseudocolor plot of the transmembrane potential V (x,y,t) showing a spiral wave at t~10s. (B) A plot of the time series of V(x,y,t) recorded from the representative point (x~2:25cm,y~2:25cm); (C) a plot versus the frequency v of the power spectrum E(v) of this time series; (D) a plot of the inter-beat interval (IBI) versus the beat number n for this time series (here we have discarded the first 10 beats to remove the initial transients. doi:10.1371/journal.pone.0045040.g013 Figure 11. Pseudocolor plots of the local membrane potential V illustrating plane-wave propagation through a transmural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. Here we consider the strong-coupling case and different régimes in the stability diagram of Fig. 10; (a.1)–(a.4) propagation in régime R1 (Pf ~0%,Ef ~{30mV ); (b.1)–(b.4) propagation in régime R2 (Pf ~5%,Ef ~{30mV); (c.1)–(c.4) propagation in régime R3 (Pf ~10%,Ef ~{30mV); (d.1)–(d.4) propagation in régime R4 (Pf ~12%,Ef ~{30mV ); (e.1)–(e.4) propagation in régime R5 (Pf ~20%,Ef ~{30mV). (For full spatiotemporal evolutions see Video S4.). doi:10.1371/journal.pone.0045040.g011 time series is a flat line, which indicates that there is no trace of activity. In contrast, the states SRSQ, MRSP, and MRSQ cannot be identified unambiguously from a quick inspection of the time series of Vm from a representative point in the simulation domain; e.g., a plot of the IBI versus n might suggest the existence of an m cycle, but only a careful analysis of the power spectrum E(v) can distinguish clearly between such an m cycle and quasiperiodic temporal evolution, with more than one, incommensurate, fundamental frequencies v1 ,v2 , etc.; furthermore, the number of spirals or rotors cannot be identified, in these cases, unless we analyze activation movies of pseudocolor plots of Vm and trace the trajectories of spiral tips (these results are in consonance with We have found both oscillatory and autorhythmic states. Although the target waves in both these cases are similar, those in the autorhythmic case have a larger amplitude than in the oscillatory case. Note, furthermore, that the states SRSP, ST, and SA can be identified merely from the time series of Vm , with data recorded from any representative point in the simulation domain: The time series for Vm , in the state SRSP, is completely periodic, so the plot of IBI versus the number n of the beat is a flat line; in the state ST this time series is obviously chaotic; in the state SA the Figure 12. The dependence of the plane-wave conduction velocity CV and wavelength l on the percentage of fibrosis Pf . Plots of (a) CV and (b) l versus Pf for the moderate-coupling case (solid black line with filled black circles), i.e., Dmf ~0:0000218 cm2/ms, and the strongcoupling case (solid blue line with unfilled blue circles), i.e., Dmf ~0:000048 cm2/ms. doi:10.1371/journal.pone.0045040.g012 PLOS ONE | www.plosone.org 13 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 14. Pseudocolor plots of the local membrane potential V illustrating spiral-wave dynamics in a mural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. We obtain six qualitatively different behaviors, namely, SRSP (Single Rotating Spiral Periodic), SRSQ (Single Rotating Spiral Quasiperiodic), MRSP (Multiple Rotating Spirals Periodic), MRSQ (Multiple Rotating Spirals Quasiperiodic), ST (Spiral Turbulence), and SA (Spiral Absorption). Illustrative pseudocolor plots of V show the time evolution of a spiral wave for (a.1)–(a.5) SRSP with Pf ~5%, (b.1)–(b.5) SRSQ with Pf ~21%, (c.1)–(c.5) MRSP with Pf ~17:8%, (d.1)–(d.5) MRSQ with Pf ~18:2%, (e.1)–(e.5) ST with Pf ~20:5%, and (f.1)–(f.5) SA with Pf ~25%. (For full spatiotemporal evolutions see Video S5.). doi:10.1371/journal.pone.0045040.g014 wave-front to the wave-back also decreases. In Fig. 18 we show, for representative values of Pf , isosurface plots of the local transmembrane potential V (x,y,t) that illustrate the time evolution of a scroll wave in three qualitatively different states, namely, SRS, MRS, and SA, which we have defined above; the spatiotemporal evolution of V (x,y,t) for these states is shown in Video S6. The states SRS and MRS have single and multiple scrolls, respectively; their temporal evolution may be periodic, quasiperiodic, or chaotic; to determine this unambiguously, we need far longer time series than we have been able to get with our computational resources. However, we can distinguish clearly between the states SRS, MRS, and SA. Given our initial conditions and the distributions of fibroblasts, the sequence of transitions in our 3D model is shown in Fig. 17(b). As we have earlier studies of spiral waves in mathematical models of cardiac tissue [43–45] without fibroblasts). Influence of diffuse fibrosis on scroll waves in 3D We consider now the dynamics of scroll waves of electrical activation in our mathematical model in the presence of fiber anisotropy and diffuse fibrosis. We start with a rectangular parallelepiped of myocytes (Pf ~0%) and the initial condition of Fig. 1 (b). We find that, even after t~4 s, the medium supports only one, temporally periodic, rotating scroll wave, which does not break up further into smaller scrolls. We call this state SRS (SingleRotating-Scroll). We now increase Pf in steps of 1% and find that, as Pf increases, this periodic, scroll-wave develops a granular texture, whose granularity increases with Pf ; the distance from the PLOS ONE | www.plosone.org 14 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 15. Time series of V illustrating high-order temporal cycles during spiral-wave propagation in a mural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. Plots of the time series of V (x,y,t), from representative points in the simulation domain, and the resulting plots of the interbeat interval IBI versus the beat number n and the power spectrum E(v) versus the frequency v illustrating temporal 2-cycles (a.1)–(a.3) for Pf ~21%, 3-cycles (b.1)–(b.3) for Pf ~21:3%), 4-cycles (c.1)–(c.3) for Pf ~16:8%), and 5-cycles (d.1)–(d.3) for Pf ~20:3%; these cycles show up most clearly in the IBI plots (a.2), (b.2), (c.2), and (d.2); but their presence can also be surmised from the time series of V (a.1), (b.1), (c.1), and (d.1) and the sharp peaks in the power spectra (a.3), (b.3), (c.3), and (d.3)). doi:10.1371/journal.pone.0045040.g015 heterogeneity, myocytes and fibroblasts. Our mathematical model introduces fibroblasts randomly, to mimic diffuse fibrosis, in the TNNP model [16,23] for human ventricular tissue; the passive fibroblasts in our model do not exhibit an action potential in the absence of coupling with myocytes; and we allow for a coupling between nearby myocytes and fibroblasts. Our in silico study is designed to explore effectively biophysically relevant ranges of the parameters that characterize myocytes, noted in the 2D case, the exact sequence in which these transitions occur depends sensitively on the initial conditions, boundary effects, and the realizations of fibroblast distributions. Discussion We have presented a comprehensive numerical study of spiraland scroll-wave dynamics in a state-of-the-art mathematical model for human ventricular tissue with fiber rotation, transmural PLOS ONE | www.plosone.org 15 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 16. Time series of V illustrating high-order temporal cycles during spiral-wave propagation in a mural slice of our 2D simulation domain with a random distribution of myocytes and fibroblasts. Plots of the time series of V(x,y,t), from the representative point (x~67:5mm,y~67:5mm), and the resulting plots of the interbeat interval IBI versus the beat number n and the power spectrum E(v) versus the frequency v illustrating temporal 6-cycles (a.1)–(a.3) for Pf ~21:1%, 7-cycles (b.1)–(b.3) for Pf ~20:7%), 9-cycles (c.1)–(c.3) for Pf ~16:6%), and 10-cycles (d.1)–(d.3) for Pf ~16:1%; these cycles show up most clearly in the IBI plots (a.2), (b.2), (c.2), and (d.2); but their presence can also be surmised from the time series of V (a.1), (b.1), (c.1), and (d.1) and the sharp peaks in the power spectra (a.3), (b.3), (c.3), and (d.3)). doi:10.1371/journal.pone.0045040.g016 fibroblasts, and their interactions. Thus, our work complements, in an important way, experimental studies of electrical-wave propagation in fibrotic cardiac tissue [4,25]; and, as we have mentioned above, it extends significantly the numerical studies initiated by Panfilov [11] and ten Tusscher, et al. [1,12,13]. Simulations by Maleckar, et al. [46] on a rabbit ventricular model suggest that the myocyte resting potential and AP PLOS ONE | www.plosone.org waveform, in the case of atrial arrhythmias, are modulated strongly by the properties and number of coupled fibroblasts, the degree of coupling, and the pacing frequency. Xie, et al. [47] have shown that a fibroblast, coupled with a myocyte, generates a gap-junction current, which flows from the myocyte to the fibroblasts and vice versa, with two main components: an early pulse of transient outward current and a 16 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions case we consider, on the random distribution of fibroblasts. However, the important qualitative points to note in our study are that (a) there is a variety of nonequilibrium states and (b) a rich sequence of transitions between them. These states can have important physical consequences. In particular, we speculate that the autorhythmic and oscillatory behaviors in the states R3 and R4 offer a possible model for ectopic foci. Thus, our studies of plane-, spiral-, and scroll-wave dynamics in our simulation domains with myocytes and fibroblasts can provide important qualitative insights into the possible effects of fibrosis on the propagation of electrical waves of activation in human ventricular tissue. In this sense, our work also builds upon the following studies: The in-vitro investigations of Miragoli, et al. [49] also suggest that fibroblasts, introduced into myocardial tissue by pressure overload or infarction, might lead to arrhythmogenesis via ectopic activity; the numerical studies of Jacquemet [50] also suggest that pacemaker-type activity can result from the coupling of cardiomyocytes with non excitable cells like fibroblasts; and Kryukov, et al. [51] have concluded, via in vitro and numerical studies of heterogeneous cardiac cell cultures and mathematical models thereof, that mixtures of excitable cells, which are initially silent, and passive cells can show transitions to states with oscillatory behavior. Interesting nonequilibrium transitions between different dynamical regimes have also been seen studied recently in a two-dimensional model for uterine tissue [52]. Our results are qualitatively in consonance with those of McDowell, et al. [22], who have used the Mahajan model [53] of the rabbit ventricular myocyte in a monodomain model in an anatomically realistic rabbit ventricular domain. In particular, they find that low densities of fibroblasts do not have a significant influence on the susceptibility to arrhythmias, moderate levels of fibroblasts increase the propensity for arrhythmias because of APD dispersion, and high fibroblast densities lead to conduction blockage. Their simulation domain is anatomically realistic whereas ours is not; however, we use the TNNP model for human cardiac tissue in contrast to the rabbit-ventricular model employed by them; furthermore, we carry out simulations at many more values of the fibroblast concentration than they do and, therefore, our simulations can uncover the details of the nonequilibrium transitions from single rotating spiral or scroll waves to the absorption state with no waves. Tanaka, et al. [54] have studied how the random distribution of fibroblasts affects the dynamics of atrial fibrillation (AF) in sheep cardiac tissue in which heart failure (HF) has been induced artificially; they have found that the number of fibrous patches is significantly larger after HF than in a control sample. They have also carried out simulation studies by using a two-dimensional human atrial model with structural and ionic remodeling that produce HF; in these simulations they demonstrate that changes in AF activation frequency and dynamics are controlled by the interaction of electrical waves with clusters of fibrotic patches. Muñoz, et al. [55] have carried out optical-mapping experiments in hetero-cellular monolayers of rat cardiac cells. Their study is designed to test whether fibroblast infiltration modifies the dynamics of spiral waves of electrical activation in such monolayers. One half of the monolayer has a randomly distributed myocyte-fibroblast mixture; the other half has a much larger concentration of myocytes (§95%) than of fibroblasts. In the former case, they find that slow (2:75 Hz), sustained re-entry is stabilized; and the wavefront propagates preferentially in the region with a high concentration of myocytes, at twice the conduction velocity (CV ) than in the region with 50% fibroblasts. Clinically, the distribution of fibroblasts, in cardiac tissue from a normal, healthy, human heart, has been found to be of the Figure 17. Representative band diagrams of states, in our 2D and in 3D studies, illustrating transitions between different spiral-wave (for 2D) and scroll-wave (for 3D) states as a function of Pf . Top panel: This band diagram shows the rich sequence of transitions, from one nonequilibrium state to another, that takes us from the state SRSP, which occurs predominantly at low values of Pf , to the state SA, which occurs at large values of Pf ; the values of Pf % are given below the band and the six states SRSP- SA are shown by a gray scale. Bottom panel: This band diagram shows the sequence of transitions, from one nonequilibrium state to another, that takes us from the state SRS, which occurs predominantly at low values of Pf , to the state SA, which occurs at predominantly large values of Pf ; the values of Pf % are given below the band and the three states SRS- SA are shown by a gray scale. The fine resolution of the transitions in 2D (top panel) cannot be achieved easily in 3D (bottom panel) without a prohibitive increase in computational costs. doi:10.1371/journal.pone.0045040.g017 later, background current during the repolarizing phase. Depending on the relative strengths of the two components, the fibroblastmyoycte coupling can alter repolarization and Cai cycling alternans, at both the cellular and tissue scales. Furthermore, in a separate study [48], they show that fibroblasts affect cardiac conduction, by creating electrotonic loading and elevating the myocyte resting potential; and they suggest that fibroblast-myocyte coupling prolongs the myocyte refractory period, which may facilitate induction of reentry in cardiac tissue with fibrosis. They have used the Luo-Rudy (I) (LRI) model and a rabbit ventricular myocyte model for their studies. Our investigation of a single MF composite, with a single myocyte coupled to Nf fibroblasts via a gap-junctional conductance Ggap , reveals five qualitatively different responses for this composite, namely, R1–R5. In R1 the response of the MF composite to an external electrical stimulus is like that of an uncoupled myocyte; in R2 this response has an additional action potential; responses R3 and R4 are autorhythmic and oscillatory, respectively; in R5 the MF composite produces a single AP after which it reaches a time-independent, non-excitable state. Our studies of 2D domains with a random distribution of fibroblasts in a myocyte background reveal that, as the percentage Pf of fibroblasts increases, the CV of a plane wave decreases, slowly at first and rapidly thereafter, until it reaches zero and there is conduction failure. If we consider spiral-wave dynamics in such a medium we find, in 2D, a variety of nonequilibrium states, temporally periodic (SRSP and MRSP), quasiperiodic (SRSQ, MRSQ), chaotic ( ST), and quiescent ( SA), and an intricate sequence of transitions between them (see Fig. 17 (a)). The analogous sequence of transitions for 3D scroll waves is given in Fig. 17 (b). As we have noted above, such transitions between nonequilibrium states in extended dynamical systems are known in a variety of problems including the onset of turbulence in pipe flow [40], dynamo transitions in magnetohydrodynamics [41], and the turbulence-induced melting of vortex crystals in two-dimensional soap films [42]. The precise sequence of such transitions often depends on initial conditions, boundary conditions, and, in the PLOS ONE | www.plosone.org 17 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Figure 18. Pseudocolor isosurface plots of the local membrane potential V illustrating scroll-wave dynamics in a mural slice of our 3D simulation domain with a random distribution of myocytes and fibroblasts. We obtain three qualitatively different behaviors, namely, SRS (Single Rotating Scroll), MRSP (Multiple Rotating Scrolls), and SA (Scroll Absorption). Illustrative pseudocolor plots with isosurface slicing of V show the time evolution of a scroll wave for (a.1)–(a.3) SRS with Pf ~1%, (b.1)–(b.3) MRS with Pf ~17%, and (c.1)–(c.5) SA with Pf ~15%. (For full spatiotemporal evolutions see Video S6.). doi:10.1371/journal.pone.0045040.g018 the gap-junction current. These cellular and ionic mechanisms may contribute to the risk of arrhythmia in fibrotic hearts. The principal limitations of our study are that we use a monodomain description for cardiac tissue and we do not use an anatomically realistic simulation domain. These lie beyond the scope of this study. However, studies by Potse, et al. [69] have compared potentials resulting from normal depolarization and repolarization in a bidomain model with those of a monodomain model; these studies show that the differences between results obtained from a monodomain model and those obtained from a bidomain model are extremely small. We intend to study our MFcomposite models in anatomically realistic domains and with their bidomain generalizations presently. A detailed study of diffuse fibrosis in an anatomically realistic rabbit ventricle is contained in Ref. [22]. In a separate study, we have also investigated [25] spiral-wave dynamics in a variant of our mathematical model that is motivated by the experiments of Refs. [3,70]. Lastly, the difference between the sizes of the myocytes and fibroblasts is accounted for, in one way, in our model, namely, by virtue of the dependence of the total cellular capacitances of these following two types: (i) long, string-type deposits of collagen or (ii) diffuse and randomly distributed patches [56,57]. With advancing age, structural remodeling occurs in the heart; this involves the proliferation of fibroblasts and the formation of interstitial collagen [56,58]. It has also been established that there is a significant correlation between increased amounts of fibrotic tissue in the heart and increased incidences of atrial and ventricular tachyarrhythmias and sudden cardiac death [59–67]. Furthermore, the partial decoupling of muscle fibers, a decrease in CV , and conduction blocks have been attributed to an increase in fibrosis [57]; and there is growing consensus that impaired electrical conduction, which can lead to the formation and breakage of spiral- and scroll- waves of electrical activity, plays an important, though perhaps not exclusive, role in arrhythmogenesis. Nguyen, et al. [68] have used a dynamic voltage-patch-clamp technique on adult rabbit ventricular myocytes, to reveal that the coupling of myocytes to myofibroblasts promotes the formation of early-after-depolarizations (EAD) as a result of a mismatch in early- versus late-repolarization reserve caused by a component of PLOS ONE | www.plosone.org 18 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Video S4 Plane-wave propagation in the 2D TNNP model in the presence of fiber anisotropy, transmural heterogeneity, randomly distributed fibroblast and strong coupling between the myocytes and the fibroblasts: We show the spatiotemporal evolution of the plane waves, via pseudocolor plots of the local transmembrane potential V (x,y,t), for (A) régime R1 (parameters as in Figs. 11 (a.1)–(a.4)), (B) régime R2 (parameters as in Figs. 11 (b.1)–(b.4)), (C) régime R3 (parameters as in Figs. 11 (c.1)–(c.4)), (D) régime R4 (parameters as in Figs. 11 (d.1)–(d.4)), and (E) régime R5 (parameters as in Figs. 11 (e.1)–(e.4)). The time interval covered is 0sƒtƒ1s, and number of frames per second is 25. (MPEG) two types of cells, because they depend on the surface areas of these cells. Aside from this, our model does not account explicitly for the differences in sizes between myocytes and fibroblasts. However, at large values of Pf , it is essential to account for fibroblast size in a more realistic way than we have. One possible way of doing this is to follow the study of Kryukov, et al. [51] in which Nf fibroblasts are allowed to couple to one myocyte; we have studied this for Nf w1 at the level of a single MF composite. The extension of this to two- and three-dimensional domains lies beyond the scope of our paper and will be taken up in a future study. Supporting Information Video S1 Plane-wave propagation in the 2D TNNP model with fiber anisotropy, randomly distributed fibroblasts, a mural section, and moderate coupling between the myocytes and the fibroblasts; panels (A), (B), (C), (D), (E), and (F), with Pf ~5%,10%,15%,20%,25%, and 30% respectively, show the spatiotemporal evolution of the plane waves in Figs. 7 (a.1), (b.1), (c.1), (d.1), (e.1), and (f.1), via pseudocolor plots of the local transmembrane potential V (x,y,t) for the time interval 0sƒtƒ0:3s, at 25 frames per second. (MPEG) Video S5 Spiral-wave dynamics in the 2D TNNP model with diffuse fibrosis. Here we show the spatiotemporal evolution of the spiral waves in Fig. 14, for the representative values of Pf considered there, via pseudocolor plots of the local transmembrane potential V (x,y,t) in the following six states: (A) a single spiral that rotates periodically SRSP, (B) a single spiral that rotates quasiperiodically SRSQ, (C) multiple spirals whose temporal evolution is periodic MRSP, (D) multiple spirals whose temporal evolution is quasiperiodic MRSQ, (E) spiral-wave turbulence ST, and (F) a state SA in which the spiral wave is absorbed at the boundaries of our simulation domain. The time interval covered is 0sƒtƒ20s, and number of frames per second is 10. (MPEG) Video S2 Plane-wave propagation, shown via pseudocolor plots of the local transmembrane potential V (x,y,t), in the 2D TNNP model with fiber anisotropy, randomly distributed fibroblasts, a mural section, and strong coupling between the myocytes and the fibroblasts for (A) régime R1 (parameters as in Figs. 8 (a.1)– (a.4)), (B) régime R2 (parameters as in Figs. 8 (b.1)– (b.4)), (C) régime R3 (parameters as in Figs. 8 (c.1)– (c.4)), (D) régime R4 (parameters as in Figs. 8 (d.1)– (d,4)), and (E) régime R5 (parameters as in Figs. 8 (e.1)– (e.5)) for the time interval 0sƒtƒ1s, at 25 frames per second. (MPEG) Video S6 Scroll-wave dynamics in the 3D TNNP model with diffuse fibrosis: We show, via isosurface plots of the local transmembrane potential V (x,y,t), the time evolution of a scroll wave in the following three states (for the representative values of Pf in Fig. 18): (A) single rotating scroll SRS, (B) multiple rotating scrolls MRS, and (C) SA, which is characterized by scroll-wave absorption at the boundaries. The time interval covered is 0sƒtƒ0:64s, and number of frames per second is 10. (MPEG) Video S3 Plane-wave propagation in the 2D TNNP model in the presence of fiber anisotropy, transmural heterogeneity, randomly distributed fibroblasts, and moderate coupling between the myocytes and the fibroblasts. We show the spatiotemporal evolution of the plane waves, via pseudocolor plots of the local transmembrane potential V (x,y,t), for (A) Pf ~5% (parameters as in Figs. 9 (b.1)), (B) Pf ~15% (parameters as in Figs. 9 (c.1)), (C) Pf ~25% (parameters as in Figs. 9 (d.1)), (D) Pf ~35% (parameters as in Figs. 9 (e.1)), and (E) Pf ~40% (parameters as in Figs. 9 (f.1)). The time interval covered is 0sƒtƒ0:3s, and number of frames per second is 25. (MPEG) 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, Microsoft Research (India), and the Robert Bosch Centre for Cyber Physical Systems (IISc) for support, and Supercomputer Education and Research Centre (SERC, IISc) for computational resources. Author Contributions Conceived and designed the experiments: RM ARN RP. Performed the experiments: RM. Analyzed the data: RM ARN RP. Contributed reagents/materials/analysis tools: RM RP. Wrote the paper: RM RP. References 5. Starobin JM, Starmer CF (1996) Boundary-layer analysis of waves propagating in an excitable medium: Medium conditions for wave-front-obstacle separation. Phys Rev E 54:430. 6. Xie F, Qu Z, Garfinkel A (1998) Dynamics of reentry around a circular obstacle in cardiac tissue. Phys Rev E 58:6355–6358. 7. Pertsov A (1997) Scale of geometric structures responsible for discontinuous propagation in myocardial tissue. In: Spooner P, Joyner RW, Jalife J. Discontinuous Conduction in the Heart. Armonk, NY: Futura Publishing Company. 8. Turner I, Huang C, Saumarez RC (2005) Numerical simulation of paced electrogram fractionation: relating clinical observations to changes in fibrosis and action potential duration. J Cardiovasc Electrophysiol 16:15161. 1. Ten Tusscher KHWJ, Panfilov AV (2007) Influence of diffuse fibrosis on wave propagation in human ventricular tissue. Europace 9, vi38vi45. 2. González-Rosa JM, Mercader N (2009) The epicardium: development, differentiation and its role during heart regeneration. Nature Reviews Cardiology (CNIC Ed.) 6:67–72. 3. Camelliti P, Borg TK, Kohl P (2005)Structural and functional characterisation of cardiac fibroblasts Cardiovascular Research 65:4051. 4. Kohl P, Camelliti P, Burton FL, Smith GL (2005) Electrical coupling of fibroblasts and myocytes: relevance for cardiac propagation. Journal of Electrocardiology 38:45–50. PLOS ONE | www.plosone.org 19 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions 38. Bondarenko VE, Rasmusson RL (2007) Simulations of propagated mouse ventricular action potentials: effects of molecular heterogeneity. Am J Physiol Heart Circ Physiol 293: H1816–H1832. 39. Bondarenko VE, Rasmusson RL (2010) Transmural heterogeneity of repolarization and Ca2+ handling in a model of mouse ventricular tissue. Am J Physiol Heart Circ Physiol 299(2): H454–H469. 40. Schneider TM, Eckhardt B, Yorke JA (2007) Turbulence transition and the edge of chaos in pipe flow. Phys Rev Lett 99:034502. 41. Sahoo G, Mitra D, Pandit R (2010) Dynamo onset as a first-order transition: Lessons from a shell model for Magnetohydrodynamics. Phys Rev E 81:036317. 42. Perlekar P, Pandit R (2010) Turbulence-induced melting of a nonequilibrium vortex crystal in a forced thin fluid film. New J Phys 12:023033. 43. Gray RA, Jalife J, Panfilov AV, Baxter WT, Cabo C, et al. (1995) Mechanisms of Cardiac Fibrillation. Science 270:1222–1223. 44. Gray RA, Pertsov AM, Jalife J (1998) Spatial and temporal organization during cardiac fibrillation. Nature 392:75–78. 45. Witkowski FX, Leon LJ, Penkoske PA, Giles WR, Spanok ML, et al. (1998) Spatiotemporal evolution of ventricular fibrillation. Nature 392:78–82. 46. Maleckar MM, Greenstein JL, Giles WR, Trayanova NA (2009) Electrotonic Coupling between Human Atrial Myocytes and Fibroblasts Alters Myocyte Excitability and Repolarization. Biophysical Journal 97(October):2179–2190. 47. Xie Y, Garfinkel A, Weiss JN, Qu Z (2009) Cardiac alternans induced by fibroblast-myocyte coupling: mechanistic insights from computational models. Am J Physiol Heart Circ Physiol 297:775–784 48. Xie Y, Garfinkel A, Camelliti P, Kohl P, Weiss JN, et al. (2009) Effects of fibroblast-myocyte coupling on cardiac conduction and vulnerability to reentry: A computational study. Heart Rhythm, 6(11), November. 49. Miragoli M, Salvarani N, Rohr S (2007) Myofibroblasts Induce Ectopic Activity in Cardiac Tissue. Circ Res 101:755–758. 50. Jacquemet V (2006) Pacemaker activity resulting from the coupling with nonexcitable cells. Phys Rev E 74:011908. 51. Kryukov AK, Petrov VS, Averyanova LS, Osipov GV, Chen W, et al. (2008) Synchronization phenomena in mixed media of passive, excitable, and oscillatory cells. CHAOS 18:037129. 52. Singh R, Xu J, Garnier NG, Pumir A, Sinha S (2012) Self-organized transition to coherent activity in disordered media. Phys Rev Lett 108:068102. 53. Mahajan A, Shiferaw Y, Sato D, Baher A, Olcese R, et al. (2008) A Rabbit Ventricular Action Potential Model Replicating Cardiac Dynamics at Rapid Heart Rates. Biophysical Journal 94, January:392–410. 54. Tanaka K, Zlochiver S, Vikstrom KL, Yamazaki M, Moreno J, et al. (2007) The Spatial Distribution of Fibrosis Governs Fibrillation Wave Dynamics in the Posterior Left Atrium During Heart Failure. Circ Res DOI: 10.1161/CIRCRESAHA.107.153858. 55. Muñoz V, Berenfeld O, Jalife J (2006) Fibroblast Infiltration Reduces Conduction Velocity and Leads to Reentry Multiplication in Cardiomyocyte Monolayers. Circulation 114:II_89. 56. De Bakker JM, Van Rijen HM (2006) Continuous and discontinuous propagation in heart muscle. J Cardiovasc Electrophysiol 17:56773. 57. Kawara T, Derksen R, De Groot JR, Coronel R, Tasseron S, et al. (2001) Activation delay after premature stimulation in chronically diseased human myocardium relates to the architecture of interstitial fibrosis. Circulation 104:306975. 58. De Bakker JM, Stein M, Van Rijen HM (2005) Three-dimensional anatomic structure as substrate for ventricular tachycardia/ventricular fibrillation. Heart Rhythm 2:7779. 59. Everett TH, Olgin JE (2007) Atrial fibrosis and the mechanisms of atrial fibrillation. Heart Rhythm 4:S247. 60. Saito T, Tamura K, Uchida D, Saito T, Togashi M, et al. (2007) Histopathological features of the resected left atrial appendage as predictors of recurrence after surgery for atrial fibrillation in valvular heart disease. Circ J 71:708. 61. Nakai T, Chandy J, Nakai K, Bellows WH, Flachsbart K, et al. (2006) Histologic assessment of right atrial appendage myocardium in patients with atrial fibrillation after coronary artery bypass graft surgery. Cardiology 108:906. 62. Strain JE, Grose RM, Factor SM, Fisher JD (1983) Results of endomyocardial biopsy in patients with spontaneous ventricular tachycardia but without apparent structural heart disease. Circulation 68:117181. 63. Segawa I, Suzuki T, Kato M, Tashiro A, Satodata R (1990) Relation between myocardial histological changes and ventricular tachycardia in cardiomyopathy: a study by 24-hour ecgmonitoring and endomyocardial biopsy. Heart Vessels Suppl 5:3740. 64. Assomull RG, Prasad SK, Lyne J, Smith G, Burman ED, et al. (2006) Cardiovascular magnetic resonance, fibrosis, and prognosis in dilated cardiomyopathy. J Am Coll Cardiol 48:197785. 65. John BT, Tamarappoo BK, Titus JL, Edwards WD, Shen WK, et al. (2004) Global remodeling of the ventricular interstitium in idiopathic myocardial fibrosis and sudden cardiac death. Heart Rhythm 1:1419. 66. Hsia HH, Marchlinski FE (2002) Characterization of the electroanatomic substrate for monomorphic ventricular tachycardia in patients with nonischemic cardiomyopathy. Pacing Clin Electrophysiol 25:111427. 67. Varnava AM, Elliott PM, Mahon N, Davies MJ, McKenna WJ (2001) Relation between myocyte disarray and outcome in hypertrophic cardiomyopathy. Am J Cardiol 88:2759. 9. Spach MS, Heidlage JF, Dolber PC, Barr RC (2007) Mechanism of origin of conduction disturbances in aging human atrial bundles: experimental and model study. Heart Rhythm 4:17585. 10. Kuijpers NH, Keldermann RH, Arts T, Hilbers P (2005) Computer simulations of successful defibrillation in decoupled and non-uniform cardiac tissue. Europace 7:16677. 11. Panfilov AV (2002) Spiral breakup in an array of coupled cells: the role of the intercellular conductance. Phys Rev Lett 88(11):118101. 12. Ten Tusscher KHWJ, Panfilov AV (2003) Influence of nonexcitable cells on spiral breakup in twodimensional and three-dimensional excitable media. Phys Rev E 68:062902. 13. Ten Tusscher KHWJ, Panfilov AV (2005) Wave Propagation in Excitable Media with Randomly Distributed Obstacles. SIAM Journal of Multiscale Modeling & Simulation 3:265–282. 14. Shajahan TK, Sinha S, Pandit R (2007) Spiral-wave dynamics depend sensitively on inhomogeneities in mathematical models of ventricular tissue. Phys Rev E 75:011929. 15. Shajahan TK, Sinha S, Pandit R (2009) The Mathematical Modelling of Inhomogeneities in Ventricular Tissue. In: Dana SK, Roy PK and J, editors. Complex Dynamics in Physiological Systems: From Heart to Brain. (Springer) pp 51–67. 16. Shajahan TK, Nayak AR, Pandit R (2009) Spiral-Wave Turbulence and its Control in the Presence of Inhomogeneities in Four Mathematical Models of Cardiac Tissue. PLoS ONE 4(3):e4738. 17. Majumder R, Nayak AR, Pandit R (2011) Scroll-Wave Dynamics in Human Cardiac Tissue: Lessons from a Mathematical Model with Inhomogeneities and Fiber Architecture. PLoS ONE 6(4): e18052. doi:10.1371/journal.pone.0018052. 18. Majumder R, Nayak AR, Pandit R (2011) An Overview of Spiral- and ScrollWave Dynamics In Mathematical Models for Cardiac Tissue. In: Tripathi ON, Ravens U, and Sanguinetti MC, editors. Heart Rate and Rhythm, (SpringerVerlag, Berlin, Heidelberg) 14:269–282; doi:10.1007/978–3-642–17575–6_14 19. Sinha S, Stein KM, Christini DJ (2002) Critical role of inhomogeneities in pacing termination of cardiac reentry. Chaos 12:893–902; doi:10.1063/ 1.1501176. 20. Takagi S, Pumir A, Pazo D, Efimov I, Nikolski V, et al. (2004) Unpinning and removal of a rotating wave in cardiac muscle. Phys Rev Lett 93: 058101. PMID: 15323732. 21. Biktashev VN, Holden AV (1998) Re-entrant waves and their elimination in a model of mammalian ventricular tissue. Chaos 8(1):48–56. 22. McDowell KS, Arevalo HJ, Maleckar MM, Trayanova NA (2011) Susceptibility to Arrhythmia in the Infarcted Heart Depends on Myofibroblast Density. Biophysical Journal 101(September):1307–1315. 23. Ten Tusscher KHWJ, Noble D, Noble PJ, Panfilov AV (2004) A model for human ventricular tissue. Am J Physiol Heart Circ Physiol 286. 24. MacCannell KA, Bazzazi H, Chilton L, Shibukawa T, Clark RB, et al. (2007) A Mathematical Model of Electrotonic Interactions between Ventricular Myocytes and Fibroblasts. Biophysical Journal 92:41214132. 25. Nayak AR, Shajahan TK, Panfilov AV, Pandit R (In press) Spiral-wave dynamics in a Mathematical Model of Human Venticular Tissue with Myocytes and Fibroblasts. unpublished manuscript. 26. Kohl P, Kamkin AG, Kiseleva IS, Noble D (1994) Mechanosensitive fibroblasts in the sino-atrial node region of rat heart: interaction with cardiomyocytes and possible role. Exp Physiol 79:943–956. 27. Rook MB, van Ginneken ACG, De Jonge B, El Aoumari A, Gros D, et al. (1992) Differences in gap junction channels between cardiac myocytes, fibroblasts, and heterologous pairs. Am J Physiol 263:C959–C977. 28. Remy-Jardin M, Remy J (2008) Integrated Cardiothoracic Imaging with MDCT, (Springer). 29. Luo CH, Rudy Y (1994) A dynamic model of the cardiac ventricular action potential. I. Simulations of ionic currents and concentration changes. Circ Res 74:1071–1096. 30. Luo CH, Rudy Y (1994) A dynamic model of the cardiac ventricular action potential. II. Afterdepolarizations, triggered activity, and potentiation. Circ Res 74:1097–1113. 31. Bernus O, Wilders R, Zemlin CW, Versschelde H, Panfilov AV (2002) A computationally efficient electro-physiological model of human ventricular cells. Am J Physiol Heart Circ Physiol 282: H2296. 32. Fenton F, Karma A (1998) Vortex dynamics in three-dimensional continuous myocardium with fiber rotation: Filament instability and fibrillation. Chaos 8(1):20–47. 33. Ferencik M, Abbara S, Hoffmann U, Cury RC, Brady TJ, et al. (2004) Left Ventricular Thin-Point Detection Using Multidetector Spiral Computed Tomography. The American Journal of Cardiology 93 April 1. 34. Qu Z, Kil J, Xie F, Garfinkel A, Weiss JN (2000) Scroll Wave Dynamics in a Three-Dimensional Cardiac Tissue Model: Roles of Restitution, Thickness, and Fiber Rotation. Biophysical Journal 78:2761–2775. 35. Stauffer D, Aharony A (1992) Introduction to Percolation Theory. (Taylor and Francis, London, 1985) 2nd ed. 36. Nishimori H (2001) Statistical Physics of Spin Glasses and Information Processing: An Introduction. Oxford University Press, Oxford, UK. 37. De Dominicis C, Giardina I (2006) Random Fields and Spin Glasses. Cambridge University Press, Cambridge, UK. PLOS ONE | www.plosone.org 20 October 2012 | Volume 7 | Issue 10 | e45040 Nonequilibrium Arrhythmic States and Transitions Propagation in the Human Heart. IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING 53(12), December. 70. Baudino TA, McFadden A, Fix C, Hastings J, Price R, et al. (2008) Cell patterning: interaction of cardiac myocytes and fibroblasts in three-dimensional culture. Microsc Microanal 14(2):117–125. 68. Nguyen TP, Xie Y, Garfinkel A, Qu Z (2011) Arrhythmogenic Consequences of Myofibroblast-Myocyte Coupling. Cardiovasc Res (2011) doi: 10.1093/cvr/ cvr292 69. Potse M, Dubé B, Richer J, Vinet A, Gulrajani RM (2006) A Comparison of Monodomain and Bidomain Reaction-Diffusion Models for Action Potential PLOS ONE | www.plosone.org 21 October 2012 | Volume 7 | Issue 10 | e45040