Table Of ContentThe effects of fibroblasts on wave dynamics in a mathematical model for human
ventricular tissue
Alok Ranjan Nayak1∗ and Rahul Pandit2†
1Robert Bosch Centre for Cyber Physical Systems,
Indian Institute of Science, Bangalore 560012, India.
2Centre for Condensed Matter Theory, Department of Physics,
Indian Institute of Science, Bangalore 560012, India.
We present systematic numerical studies of electrical-wave propagation in two-dimensional (2D)
andthree-dimensional(3D)mathematicalmodels, forhuman, ventriculartissuewithmyocytecells
that are attached (a) regularly and (b) randomly to distributed fibroblasts. In both these cases
we show that there is a parameter regime in which single rotating spiral- and scroll-wave states
(RS) retain their integrity and do not evolve to a state ST that displays spatiotemporal chaos
6 and turbulence. However, in another range of parameters, we observe a transition from ST to
1 RS states in both 2D or 3D domains and for both cases (a) and (b). Our studies show that the
0 ST-RS transition and rotation period of a spiral or scroll wave in the RS state depends on (i) the
2
coupling strength between myocytes and fibroblasts and (ii) the number of fibroblasts attached to
n myocytes. Weconcludethatmyocyte-fibroblastcouplingstrengthandthenumberoffibroblastsare
a more important for the ST-RS transition than the precise way in which fibroblasts are distributed
J over myocyte tissue.
1
1 PACSnumbers: 87.19.Hh,89.75.Kd,05.45.-a
Keywords: MathematicalModel,Fibroblasts,Action-potential-duration-restitution,Spiral-wavedynamics
]
O
T I. INTRODUCTION 16]. Zlochiver, et al. [15] have shown the expression
of Cx43 between fibroblasts and myocytes in a mono-
.
o Approximately 16% of all deaths in the industrialized layer of myocytes and fibroblasts of neonatal rats; Mi-
i
b world are caused by cardiac arrhythmias like ventricular ragoli, et al. [16] have reported that Cx43 and Cx45
- tachycardia (VT) and ventricular fibrillation (VF) [1, 2]. are expressed among fibroblasts and between fibrob-
q
[ There is a broad consensus that the analogs of VT and lasts and myocytes in cultured fibroblasts coated over
VF in mathematical models for cardiac tissue are, re- rat-ventricular-myocyte strands. Both experimental and
1 spectively, (a) a single rotating spiral or scroll wave of computationalstudieshaveshownthatsuchcouplingbe-
v electrical activation and (b) spiral-wave or scroll-wave tweenmyocytesandfibroblastsenhanceselectrical-signal
7
turbulence, which displays broken electrical waves and propagation in cardiac tissue [7, 15–18]; this enhance-
9
3 spatiotemporalchaos[3–6]. Thus,itisveryimportantof ment increases with Nf, the number of fibroblasts that
2 study such spiral and scroll waves to develop an under- are attached to a myocyte. In cell-culture experiments,
0 standing of life-threatening cardiac arrhythmias. Such Miragoli, et al. [16] have found that the conduction ve-
1. studies are truly interdisciplinary because they require locity (CV) decreases with an increase in the density of
0 inputs from biology, bio-medical engineering, cardiology, fibroblasts in cultured strands of neonatal-rat myocytes
6 on the one hand, and physics, nonlinear dynamics, and coated by fibroblasts; studies by McSpadden, et al. [17]
1 numerical methods, on the other. We use methods from have found that CV decreases as the fibroblast number
v: these areas to solve the nonlinear, partial-differential increases on the top of a myocyte layer in a monolayer
i equations for state-of-the-art mathematical models for of neonatal-rat cardiac myocytes, which are electrotoni-
X cardiac tissue with myocytes and fibroblasts. We build callyloadedwithalayerofcardiacfibroblasts. Zlochiver,
r on our earlier studies of this problem [7, 8] to elucidate et al. [15] have shown that CV decreases as (i) the gap-
a
the role of two different forms of myocyte-fibroblast cou- junctional conductance increases or (ii) the fibroblasts
plingsonspiral-andscroll-wavedynamicsinsuchmodels density increases in their experiments with fibroblasts of
by using theoretical ideas from spatiotemporal chaos. neonatal rats; they have also obtained similar result in
It is useful to begin with an overview of experimen- their computational studies in a two-dimensional (2D)
tal and computational studies of cardiac myocytes and sheet of myocyte tissue in the dynamic Luo-Rudy (LRd)
fibroblasts. Cardiac fibroblasts, which are inexcitable model [19, 20] by inserting fibroblasts. Computational
cells, often multiply and connect with cardiac myocytes studies by Xie, et al., [18] have shown that CV decreases
during fibrosis [9–11], a process of cardiac-tissue heal- as they increase the gap-junctional coupling or the fi-
ing after a myocardial infarction. Fibroblasts in cell cul- broblasts density in a 2D LR1 myocyte model [21], with
ture and in intact tissue can couple with myocytes by either randomly attached or randomly inserted fibrob-
expressing either the connexins-43 (Cx43) or Cx45 [12– lasts. In simulations with the 2D TNNP04 model (due
to ten Tusscher, et al. [22]), Nayak, et al. [7] have found
that CV either decreases or increases, with attached fi-
broblasts, as they increase the gap-junctional coupling.
∗ [email protected] Theexperimentalandcomputationalinvestigationsmen-
† [email protected]
2
tionedaboveshowthatboththegap-junctionalcoupling hereC isthetotalcellularcapacitanceofamyocyte,V
m m
and N enhance CV and, therefore, they can play a cru- isthemyocytetransmembranepotential,i.e.,thevoltage
f
cial role in spiral- and scroll-wave dynamics in mathe- difference between intra- and extra-cellular spaces, and
matical models for cardiac tissue. I is the sum of all the ionic currents that cross the my-
m
We develop and investigate two models with differ- ocytecellmembrane;C ,V ,andI are,respectively,the
f f f
ent arrangements of fibroblasts that are attached to my- total cellular capacitance, the transmembrane potential,
ocytes. In the first arrangement there is a regular, spa- andthesumofallioniccurrentsforthefibroblast;N (x)
f
tially periodic attachment of fibroblast, whereas, in the is the number of identical fibroblasts attached to a my-
second arrangement, fibroblasts are attached randomly ocyteinoursimulationdomainatthepointx;andI ,G ,
j j
tomyocytes. Ourstudyhasbeendesignedtounderstand and D are, respectively, the gap-junctional current, the
the effects of fibroblast organization, fibroblast density, myocyte-fibroblastsgap-junctionalconductance, andthe
andthemyocyte-fibroblastcouplingonspiral-andscroll- diffusion coefficient that is related to the gap-junctional
wavedynamics. Weusetwoparametersetsformyocytes. conductance between myocytes.
The first set leads to a stable rotating spiral or scroll For myocytes, we use the state-of-the-art mathemati-
(RS) wave; the second leads to spiral- or scroll-wave- cal model for human ventricular tissue developed by ten
turbulence (ST) states in an isolated myocyte domain. Tusscher and Panfilov (the TP06 model) [25]. In the
By investigating an ST state in the presence of fibrob- TP06 model the total ionic current is
lasts, we observe that both models, with regularly and
randomly attached fibroblasts, show transitions from an Im=INa+ICaL+Ito+IKs+IKr+IK1 (4)
ST to an RS state, depending on the myocyte-fibroblast +I +I +I +I +I +I ,
NaCa NaK pCa pK bNa bCa
coupling G and the maximum number N of fibroblasts
j f where I is the fast, inward Na+ current, I the L-
attached to a myocyte in our simulation domain. We Na CaL
type, slow, inward Ca2+ current, I the transient, out-
findthat, onceSTisconvertedtoRS,thespiralorscroll to
wardcurrent,I theslow,delayed,rectifiercurrent,I
rotationperiodincreasesasweincreaseG andN . Our Ks Kr
j f
the rapid, delayed, rectifier current, I the inward, rec-
studywithanRSstateandfibroblastsshowsthatanRS K1
tifier K+ current, I the Na+/Ca2+ exchanger cur-
remains unchanged in both models with regularly and NaCa
rent,I theNa+/K+ pumpcurrent,I andI the
randomly attached fibroblast; and the rotation period NaK pCa pK
plateau Ca2+ and K+ currents, and I and I the
increases as we increase G and N . bNa bCa
j f background Na+ and Ca2+ currents, respectively. The
The remainder of this paper is organized as follows.
full sets of equations for this model, including the ODEs
Section II is devoted to a description of our model and
for the ion-channel gating variables and the ion dynam-
the numerical methods we use. Section III is devoted
ics, are given in Refs. [7, 26].
to our results. Section IV contains a discussion of the
We follow MacCannell, et al. [27] to model the fibrob-
significance of our results.
lasts as passive elements. The fibroblast ionic current I
f
is
II. MODEL AND METHODS
I =G (V −E ), (5)
f f f f
InthisSection,wedescribethedetailsofourmyocyte- where G and E are, respectively, the conductance and
f f
fibroblast models for two-dimensional (2D) and three- the resting membrane potential for the fibroblast.
dimensional (3D) tissue. We also explain the numerical- Physical units in our model are as follows: time t is
simulation techniques that we use to solve the partial in milliseconds (ms), V and V are in millivolts (mV),
m f
differential equations (PDEs) that comprise our mathe- C and C are in picofarads (pF), I and I are in
m f m f
maticalmodels. Wealsodiscussthemethodsthatweuse picoamperes (pA), G is in nanoSiemens (nS), E is in
f f
to analyze the data from our numerical simulations. mV, G is in nS, and D is in cm2/ms.
j
We study models with (a) regularly attached fibrob-
lasts and (b) randomly attached fibroblasts. In case (a),
A. Model Nf(x) = Nf for all site x in our simulation domain. In
case(b)wechooseN (x)randomlyateachsitex;N (x)
f f
can be any integer from 0 to N , with equal probability
The 2D and 3D myocyte domains, with attached fi- f
for any one of these values.
broblasts, can be modeled by the following PDEs and
ordinary-differential-equations (ODEs) [23, 24]:
∂V −I +N (x)I B. Methods
m = m f j +D∇2V , (1)
m
∂t C
m
∂V −I −I Our 2D and 3D simulation domains are, respectively,
f f j
= , (2) squares (1024×1024 grid points) and rectangular paral-
∂t C
f lelepipeds (1024×1024×8 grid points). We use 5-point
where and 7-point stencils for the Laplacian in 2D and 3D, re-
spectively, and a finite-difference scheme with step sizes
I =G (V −V ); (3) δx=δy=0.25 mm in 2D, and δx=δy=δz =0.25 mm in
j j f m
3
3D, i.e., our simulation domains are 256×256 mm3 (in conductance,is0.00219nS/pF;and(e)τ ,thetimecon-
f
2D) and 256×256×2 mm2 (in 3D). For time marching stant of the f gating variable that is associated with the
we use a forward-Euler scheme with δt = 0.02 ms. We I current, is increased 2 times compared to its value
CaL
use Neumann (no-flux) conditions at the boundaries of intheTP06model[7,25,26]. Ourfibroblastsparameters
our simulation domain. are as follows: C =6.3 pF, G =4 nS, E =−49.0 mV,
f f f
Fornumericalefficiency,wehavecarriedoutoursimu- and Gj in the range 0≤Gj ≤6 nS [7, 30].
lations on parallel computers, with an MPI code that we To obtain spiral and scroll waves we use the S1-S2
have developed for the TP06 model. Our code divides cross-fieldprotocol[26,31,32]. Weapplyastimulus(S1)
the 2D (or 3D) simulation domain into n columns (or of strength 150 pA for 3 ms to the left boundaries of our
slabs)alongthex-directionofthedomain, i.e., eachpro- simulationdomains,toformaplanewave. Wethenapply
cessor carries out the computations for (1024/n)×1024 thesecond(S2)stimulus,withthesamestrengthanddu-
and (1024/n)×1024×16 grid points, respectively, for 2D rationastheS1stimulus,fromthebottomboundaryand
and 3D domains. To compute the Laplacian at the in- with 0 mm≤y≤125 mm in 2D, and 0 mm≤y≤125 mm
terface of processor boundaries, we use two extra grid and 0 mm≤z ≤2 mm in 3D. This protocol leads to the
lines (or surfaces), which can send and receive the data formation of spiral and scroll waves, respectively, in our
from left- and right-neighbor processors. The Neumann 2D and 3D domains.
boundary condition is taken care of by adding an extra To examine the spatiotemporal evolution of our sys-
layer of grid points on the boundaries of the simulation tem, we obtain pseudocolor or isosurface plots of Vm,
domain of each processor. time series of Vm, from representative points (x = 125
Reference[28]suggeststhatwemusthaveDδt/(δx2)< mm, y = 125 mm for 2D, and x = 125 mm, y = 125
1/2d for numerical stability, where d is the dimension of mm, z=1.25 mm for 3D), which we mark with an aster-
the simulation domain. For the TP06 model, with dif- isk (∗) in all pseudocolor plots of Vm. We examine the
fusion coefficients D = 0.00154 cm2/ms [25], time step inter-beat interval (IBI), by using this time series with
δt = 0.02 ms, and space step δx = 0.25 mm, the value 4.4×105 data points; the IBI is the interval between two
of Dδt/(δx)2 is ≃ 0.05; for the TP06 model, the quan- successivespikesinthistimeseries. Weobtainthepower
tity 1/2d = 0.25 and ≃ 0.17, for the 2D and 3D do- spectra E(ω), of the time series of Vm, by using 2×105
mains, respectively, i.e., we have numerical stability be- datapoints; toeliminatetransientsweremovetheinitial
cause Dδt/(δx)2<1/2d. 2.4×105 data points. To obtain the rotation period T
of a spiral, in an RS state, we average over the last 5
We check the accuracy of our numerical scheme, as in
rotations of that RS.
Ref. [28], by varying both δt and δx in a cable-type do-
mainofmyocytes[7,26]andbymeasuringCVofaplane
wave, which is injected into the domain by stimulating
its left boundary for 3 ms with a stimulus of strength III. RESULTS
150pA.Withδx=0.25mm, CVincreasesby1.1%when
wechangeδtfrom0.02to0.01ms;ifwedecreaseδxfrom In subsection IIIA, we begin by studying spiral-wave
0.25to0.15mm,withδt=0.02ms,thenCVincreasesby dynamics in a 2D domain of myocytes without fibrob-
3.1%; these changes are similar to those found in other lasts; we then introduce fibroblasts, either regularly or
studies [22, 26, 28, 29]. randomly, and examine the effects they have on spiral-
Although the numerical method we use satisfies both wave dynamics. Subsection IIIB contains the results of
numerical-stability and accuracy conditions, an inappro- our studies of scroll-wave dynamics in our 3D simulation
priately large δx can give irregular wavefront-curvature, domain.
as a consequence of numerical artifacts [7, 26, 28]; this
leads to unphysical wave dynamics. We check that our
results are free from such numerical artifacts by inves- A. Spiral-wave dynamics in our 2D model
tigating the spatiotemporal evolution of an expanding
wave front that emerges from a point stimulus. We find In Fig. 1(A), we show a pseudocolor plot of V , at
m
that fronts of the expanding wave do not deviate sub- time t = 8.8 s, for the parameter set P1, in our 2D sim-
stantially from circles, when we apply a point stimulus ulation domain without fibroblasts; the initial condition
of strength 450 pA for 3 ms at the center of the domain. evolves to a state with a single rotating spiral (RS). The
We use two parameter sets P1 and P2 for myocytes local time series of V (x,y,t), from the representative
m
to obtain, respectively, a stable rotating spiral (RS) and point shown by the asterisk in Fig. 1(A)), is given in
a spiral-turbulence (ST) states in our 2D simulation do- Fig. 1(B) for 0 s ≤ t ≤ 8.8 s; a plot of the IBI versus
main, and a stable rotating scroll or scroll-wave turbu- the beat number is given in Fig. 1(C), which shows that,
lence in our 3D domain. The parameter set P1 is the after initial transients, the spiral wave rotates periodi-
original one used in the TP06 model [7, 25, 26]. In the cally with an average rotation period T ≃ 212 ms. The
P2 parameter set, we use the following parameters, with power spectrum E(ω) in Fig. 1(D) has discrete peaks at
all other parameters the same as in the original TP06 thefundamentalfrequencyω ≃4.75Hzanditsharmon-
f
model: (a) G , the I conductance, is 0.172 nS/pF; ics. The periodic time series of V , the flattening of the
Kr Kr m
(b)G ,theI conductance,is0.441nS/pF;(c)G , IBI,andthediscretepeaksinE(ω)demonstratethatthe
Ks Ks pCa
theI conductance,is0.8666nS/pF;(d)G ,theI time evolution of the spiral wave, with the P1 parame-
pCa pK pK
4
FIG.1. Therotating-spiral(RS)andspiral-turbulence(ST)
states in a 2D domain in the absence of fibroblasts. For the
parameter set P1, the pseudocolor plot of V in (A), the
m
periodic nature of the local time series for V from the rep-
m
resentativepoint(markedbyanasterisk∗in(A))in(B),the
flattening IBI with an average rotation period T ≃212 ms in
(C), and the discrete peaks in the power spectrum with the
fundamentalfrequencyωf =4.75Hzanditsharmonicsin(D) FIG.2. VariousRSandSTstatesinour2Ddomainforthe
characterize the RS state. The exact analogs of plots (A)- regularlyattachedfibroblastmodel. ((A1)-(B4))Pseudocolor
(D) are shown, respectively, in (E)-(H) for the P2 parameter plots of V , with G =4 nS for N =1,2,4,and 6, illustrate
m j f
set; the irregular local time series, the fluctuating behavior thattheRSstate(parametersetP1)remainsinanRSstateas
of the IBI and the broad-band nature of the power spectrum N increases(firstrow). However,anSTstate(parameterset
f
characterize the ST state. P2)showsatransitiontoanRSstateasN increases(second
f
row). ((C)-(D))TherotationperiodTofRSincreasesasN
f
increases,forafixedvalueofG ,andvice-versa,forboththe
j
P1 and P2 parameter sets.
ter set, is periodic. In Figs. 1(E)-(H), we show the exact
analogs of Figs. 1(A)-(B) for the P2 parameter set. The
non-periodic local time series in Fig. 1(F), the fluctuat-
ingIBIinFig.1(G),andthebroad-bandnatureofE(ω)
inFig.1(H)arecharacteristicoftheSTstate. Thepseu-
slope of a myocyte-fibroblast composite at the cellular
docolorplotinFig.1(E),attimet=8.8s,showssuchan
level [8, 33]. Once the ST state is suppressed, a sin-
STstate, whicharisesfromthesteepslopeoftheaction-
gle spiral in an RS state rotates periodically as shown
potential-duration-restitution (APDR) plot [8, 25]. In
Figs. 2((B2)-(B4)). In Figs. 2(C) and (D), we plot, re-
summary, then, in the absence of fibroblasts, the param-
spectively, the rotation period T of a spiral wave in an
eter sets P1 and P2 lead, respectively, to (a) an RS state
RS state versus N , for different values of G , for the
f j
and (b) an ST state in our 2D simulation domain.
P1 and P2 parameter sets. We find that T increases as
We now examine the effects of fibroblasts on spiral-
(i) N increases, with a fixed value of G , and (ii) G
f j j
wave dynamics in both RS and ST states. We begin
increases, with a fixed value of N . This increase of T
f
our investigation with the P1 parameter set and with
is a consequence of the decrease of CV that is associated
the regularly attached fibroblast model with 1 ≤ N ≤ 6
f with a decrease of the upstroke velocity of a myocyte-
and 1 nS ≤ G ≤ 8 nS; the remaining parameters for
j fibroblast composite AP in its depolarization phase, as
fibroblasts are as in subsection IIB.
shown in Refs. [7, 8, 18]. Furthermore, we observe that
In Figs. 2(A1)-(A4) we show pseudocolor plots of V ,
m the minimum value of N , required for the ST-RS tran-
f
at time t=8.8 s, with the P1 parameter set for regularly
sition, decreases as G increases, for the P2 parameter
j
attached fibroblast model with G = 4 nS and different
j set.
valuesofN . TheanalogsofFigs.2(A1)-(A4)areshown
f
in Figs. 2(B1)-(B4) for the P2 parameter set. The plots We focus next on spiral-wave dynamics, with P1 and
inthefirstrowofFig.2showthattheRSstate,whichwe P2 parameter sets, in our randomly attached fibroblast
obtain in the absence of fibroblasts, does not evolve into model. In Fig. 3, we show the exact analogs of Fig. 2,
anSTstate; however,thespiral-armwidthW decreases butnowfortherandomlyattachedfibroblastmodel. The
d
as we increase N , for a fixed value of G (first row of pseudocolorplotsinFigs.3(A1)-(A4)showthattheran-
f j
Fig. 2). We define W to be the difference of the radial domness in attaching fibroblasts does not lead to an RS-
d
distance between the wave front and the wave back of STtransitionfortheP1parameterset. However, inspite
a spiral arm, whose center is located at the spiral core. of the randomness in the arrangement of fibroblasts, we
Such a decrease of W is related to the shortening of the observe an ST-RS transition for the P2 parameter set
d
action-potential-duration (APD) of a myocyte-fibroblast (see Figs. 3(B1)-(B4)), which is qualitatively similar to
composite that has been discussed in Refs. [7, 8]. theST-RStransitionintheP1case(comparethesecond
FortheP2parameterset,weobserveatransitionfrom rows of Figs. 2 and 3). However, the minimum value of
an ST to an RS state as we increase N for a fixed value N , required for an ST-RS transition, is higher for the
f f
of G (second row of Fig. 2). Such an ST-RS transition randomlyattachedfibroblastmodelthanintheregularly
j
istheconsequenceofthesuppressionofthesteepAPDR attached case (compare Figs. 2(D) and 3(D)).
5
Parameter set P1 Parameter set P2
350 300
s) 300 (C) GGjj==14 nnSSs) 250 (D)
(m 250 (m
200
T 200 T
150 150
0 2 4 6 0 2 4 6
N N
f f
FIG.3. VariousRSandSTstatesinour2Ddomainforthe
randomlyattachedfibroblastmodel(thisfigureistheanalog FIG. 4. The rotating-scroll and scroll-wave-turbulence
of Fig. 2). The results are qualitatively similar to those in states, in our 3D simulation domain of size 256 × 256 ×
Fig. 2 for the regularly attached fibroblast case. Note that 2mm3,fortheregularlyattachedfibroblastmodel,areshown
the minimum value of Nf, for a fixed value of Gj, for the in ((A1)-(B4)) via isosurface plots of Vm. The myocyte-
ST-RS transition, is higher compared to that in Fig. 2(D). fibroblast coupling strength Gj =4 nS. The scroll-arm width
of a rotating scroll, with the P1 parameter set, decreases as
N increases (first row). The scroll-wave turbulence, associ-
f
ated with the P2 parameter set, is converted to a rotating
scroll as N increases (second row). ((C)-(D)) Plots of the
f
B. Scroll-wave dynamics in our 3D model
rotation period T of a scroll wave in a rotating-scroll state,
for the P1 and P2 parameter sets; for both parameter sets T
We turn now to a systematic study of scroll-wave dy- increasesasGj increases;notethat,fortheST-RStransition,
namics in our 3D simulation domain. For both the P1 the minimum value of Nf decreases as Gj increases.
and P2 parameter sets and both regularly and randomly
attached fibroblast models, we carry out simulations to
studythedependenceofscroll-wavedynamicsonN and
f
G . We present our numerical results below.
j
In Figs. 4(A1), (A2), (A3) and (A4), we show, respec-
tively, isosurface plots of V , at time t = 8.8 ms, for
m
theP1parametersetinourregularlyattachedfibroblast
model with G = 4 nS and N = 0 (i.e., isolated my-
j f
ocytes), N = 1, N = 2, and N = 4. In the absence
f f f
of fibroblasts, i.e., N = 0, the P1 parameter set dis-
f
plays a rotating scroll wave with fundamental frequency
ω ≃5 Hz and rotation period T ≃ 201 ms; this is consis-
f
tent,becauseω ≃1/T. InFig.4(C),weplotTversusN
f f
for G = 1 nS (●) and 4 nS (▴). We find that T increases
j
as we increase (i) N , for a fixed value of G , or (ii) G , Parameter set P1 Parameter set P2
f j j 215 190
for a fixed value of N . In Figs. 4((B1)-(B4)) and (D),
f (C) Gj=1 nS (D)
we show, respectively, the exact analogs of Figs. 4((A1)- s) 210 Gj=4 nSs)
(A4)) and (C), for the P2 parameter set. In the absence (m (m 180
205
of fibroblasts and for the P2 parameter set, we obtain T T
a scroll-wave-turbulence state (Fig. 4(B1)); this scroll- 200 170
wave turbulence is converted to a rotating scroll if we 0 2 N 4 6 0 2 N 4 6
f f
have N >1 (second row of Fig. 4). Once the scroll-wave
f
turbulence state is suppressed, a rotating scoll rotates
FIG. 5. The rotating-scroll and scroll-wave-turbulence
with a period T, which increases as we increase N for
f states,inour3Dsimulationdomainfortherandomlyattached
a fixed value of G , and vice-versa (Fig. 4(D)). Further-
j fibroblast model; the exact analog of Fig. 4. The results are
more, from Fig. 4(D), we find that the minimum value qualitativelysimilartothosefortheregularlyattachedfibrob-
of Nf, required for the ST-RS transition, is 4 and 2, re- lastcase. NotethattheST-RStransitionNf value,forafixed
spectively,forG =1nS(●)and4nS(▴). Theisosurface value of G , is higher than its counterpart in Fig. 4(D).
j j
plots in Fig. 4 show that the width W of a scroll-wave
d
arm in the rotating-scroll state decreases as we increase
N for both the P1 and P2 parameter sets. The mecha-
f
6
nisms of the ST-RS transition, and increase of T and a spiral-wavebreakupoccurs,inanLR1model,becauseof
decrease of W , as we increase N and G , are the same randomly diffuse fibroblasts in a localized area of a sim-
d f j
as those we have found in our 2D studies. ulation domain; Zlochiver, et al. [15] have shown from
InFig.5weshowtheexactanalogofFig.4fortheran- their experiments and simulations that a rotating spiral
domly attached fibroblast model, with both the P1 and becomes unstable and, finally, spiral breakup occurs, as
P2 parameter sets. Our scroll-wave results here are sim- they increase the percentage of diffuse fibroblasts. Ma-
ilar to those for the case of regularly attached fibroblast jumder, et al. [34] have shown from their numerical ex-
model. From Fig. 5(D), we find that the minimum value periments that a transition from an RS to various ST
of N , for the ST-RS transition, is 5 and 3, respectively, states occurs depending on the percentage of fibroblasts
f
for G =1 nS (●) and 4 nS (▴). Note that this minimum in their simulation domain. In our attached-fibroblast
j
value of N is higher for the randomly attached fibrob- model studies, fibroblasts do not inhibit wave propa-
f
last model than it is for the regularly attached fibroblast gation [7]; however, fibroblasts attached to a myocyte
model (compare Figs.4(D) and 5(D)). can lower the steepness of the APDR curve, depending
on the values of G and N [8, 33]. Such a lowering
j f
of the steep slope of the APDR eliminates spiral- and
IV. CONCLUSIONS scroll-wave turbulence in our 2D and 3D simulation do-
mains [35–38]. Therefore, we observe an ST-RS transi-
Wehavepresentedthemostextensivenumericalstudy tion. Earlier studies in Ref. [33] have observed ST-RS
carried out so far of the effects of fibroblasts on spiral- spiral-wave transitions because of a suppression of the
and scroll-wave dynamics in a mathematical model for steep portion of the APDR slope in a 3D model consist-
humanventriculartissuewithfibroblasts,attachedregu- ing of myocytes, fibroblasts, and extracellular space by
larlyorrandomlytomyocytes. Ournumericalstudyhas using the LR1 model [21]. However, those studies have
been designed to uncover the role of (i) the organization not investigated the spiral- and scroll-wave transition as
of fibroblasts in ventricular tissue ( i.e., to compare reg- a function of Nf and Gj. Our study shows that both
ular and random arrangements), (ii) myocyte-fibroblast Nf and Gj are important factors during the fibrosis pro-
coupling G , and (iii) the density of fibroblasts, i.e., the cess [10, 39].
j
maximum number of fibroblasts N attached to a my- We suggest that our results from in silico studies can
f
ocyte. One of the principal results of our studies is that be verified in in vitro experiments. Furthermore, by us-
spiral- and scroll-wave dynamics depend only slightly on ing advanced cell-culture techniques [40–42], our 2D and
the details of the organization of fibroblasts in ventricu- 3D numerical results can be tested easily in cell-culture
lar tissue. However, the ST-RS transition, the stability experiments.
of spiral- and scroll-wave turbulence, the rotation period
ofarotatingspiralandscroll,andthewidthofarotating
spiral and scroll arms, depend sensitively on N and G .
f j
ACKNOWLEDGMENTS
Earlier studies have investigated the effects of fibrob-
lasts on spiral-wave dynamics by introducing randomly
diffuse fibroblasts in a myocyte domain [15, 34]. Such We thank the Department of Science and Technology
randomlydiffusefibroblastsinamyocytedomaininhibit (DST),India,theUniversityGrantsCommission(UGC),
electrical-wave propagation, and initiate spiral-wave tur- India, and the Robert Bosch Centre for Cyber Physical
bulence state. Studies by Xie, et al [18] have found that Systems (RBCCPS), IISc, for support.
[1] A. S. Go, D. Mozaffarian, V. L. Roger, E. J. Benjamin, [9] G. Krenning, E. M. Zeisberg, and R. Kalluri, Journal of
J. D. Berry, M. J. Blaha, S. Dai, E. S. Ford, C. S. Fox, cellular physiology 225, 631 (2010).
S. Franco, et al., Circulation 129, e28 (2014). [10] A.BiernackaandN.G.Frangogiannis,Aginganddisease
[2] D. P. Zipes and H. J. Wellens, Circulation 98, 2334 2, 158 (2011).
(1998). [11] C. A. Souders, S. L. Bowers, and T. A. Baudino, Circu-
[3] J.M.Davidenko,A.V.Pertsov,R.Salomonsz,W.Bax- lation research 105, 1164 (2009).
ter, and J. Jalife (1992). [12] M. Rook, A. Van Ginneken, B. e. de Jonge,
[4] A. M. Pertsov, J. M. Davidenko, R. Salomonsz, W. T. A. El Aoumari, D. Gros, and H. Jongsma, American
Baxter,andJ.Jalife,Circulationresearch72,631(1993). JournalofPhysiology-CellPhysiology263,C959(1992).
[5] R.A.Gray,A.M.Pertsov,andJ.Jalife,Nature392,75 [13] P.Camelliti,G.P.Devlin,K.G.Matthews,P.Kohl,and
(1998). C. R. Green, Cardiovascular research 62, 415 (2004).
[6] J.Jalife,R.A.Gray,G.E.Morley,andJ.M.Davidenko, [14] L. Chilton, W. R. Giles, and G. L. Smith, The Journal
Chaos: AnInterdisciplinaryJournalofNonlinearScience of physiology 583, 225 (2007).
8, 79 (1998). [15] S. Zlochiver, V. Munoz, K. L. Vikstrom, S. M. Taffet,
[7] A. R. Nayak, T. Shajahan, A. Panfilov, and R. Pandit, O. Berenfeld, and J. Jalife, Biophysical journal 95, 4469
PLoS One 8, e72950 (2013). (2008).
[8] A. R. Nayak and R. Pandit, Physical Review E 92, [16] M. Miragoli, G. Gaudesius, and S. Rohr, Circulation re-
032720 (2015). search 98, 801 (2006).
7
[17] L.C.McSpadden, R.D.Kirkton, andN.Bursac, Amer- [31] A. R. Nayak and R. Pandit, Frontiers in physiology 5,
ican Journal of Physiology-Cell Physiology 297, C339 207 (2014).
(2009). [32] R.Majumder,A.R.Nayak,andR.Pandit,inHeartRate
[18] Y. Xie, A. Garfinkel, P. Camelliti, P. Kohl, J. N. Weiss, and Rhythm (Springer, 2011), pp. 269–282.
and Z. Qu, Heart Rhythm 6, 1641 (2009). [33] V. S. Petrov, G. V. Osipov, and J. Kurths, Chaos: An
[19] C.-h. Luo and Y. Rudy, Circulation research 74, 1071 InterdisciplinaryJournalofNonlinearScience20,045103
(1994). (2010).
[20] C.-H. Luo and Y. Rudy, Circulation Research 74, 1097 [34] R. Majumder, A. R. Nayak, and R. Pandit, PloS one 7,
(1994). e45040 (2012).
[21] C.-h. Luo and Y. Rudy, Circulation research 68, 1501 [35] Z.Qu,J.N.Weiss,andA.Garfinkel,PhysicalReviewE
(1991). 61, 727 (2000).
[22] K. Ten Tusscher, D. Noble, P. Noble, and A. Panfilov, [36] Z. Qu, F. Xie, A. Garfinkel, and J. N. Weiss, Annals of
American Journal of Physiology-Heart and Circulatory biomedical engineering 28, 755 (2000).
Physiology 286, H1573 (2004). [37] J.N.Weiss,P.-S.Chen,Z.Qu,H.S.Karagueuzian,and
[23] J. Keener and J. Sneyd, Mathematical physiology, inter- A. Garfinkel, Circulation Research 87, 1103 (2000).
disciplinary applied mathematics 8 (1998). [38] A. Garfinkel, Y.-H. Kim, O. Voroshilovsky, Z. Qu, J. R.
[24] A.V.Panfilov,A.V.Holden,etal.,Computational biol- Kil, M.-H. Lee, H. S. Karagueuzian, J. N. Weiss, and
ogy of the heart (Wiley, 1997). P.-S.Chen,ProceedingsoftheNationalAcademyofSci-
[25] K. Ten Tusscher and A. Panfilov, Am J Physiol Heart ences 97, 6061 (2000).
Circ Physiol 291, H1088 (2006). [39] T.P.Nguyen,Z.Qu,andJ.N.Weiss,Journalofmolec-
[26] A. R. Nayak, Ph.D. thesis, Indian Institute of Science ular and cellular cardiology 70, 83 (2014).
(2013). [40] Y. Haraguchi, T. Shimizu, T. Sasagawa, H. Sekine,
[27] K.A.MacCannell,H.Bazzazi,L.Chilton,Y.Shibukawa, K.Sakaguchi,T.Kikuchi,W.Sekine,S.Sekiya,M.Yam-
R. B. Clark, and W. R. Giles, Biophysical journal 92, ato, M. Umezu, et al., Nature protocols 7, 850 (2012).
4121 (2007). [41] T.Shimizu,M.Yamato,Y.Isoi,T.Akutsu,T.Setomaru,
[28] R. Clayton and A. Panfilov, Progress in biophysics and K. Abe, A. Kikuchi, M. Umezu, and T. Okano, Circula-
molecular biology 96, 19 (2008). tion research 90, e40 (2002).
[29] T. Shajahan, A. R. Nayak, and R. Pandit, PLoS One 4, [42] T. A. Baudino, A. McFadden, C. Fix, J. Hastings,
e4738 (2009). R. Price, and T. K. Borg, Microscopy and microanalysis
[30] Y. Xie, A. Garfinkel, J. N. Weiss, and Z. Qu, American 14, 117 (2008).
Journal of Physiology-Heart and Circulatory Physiology
297, H775 (2009).
