Spiral-wave dynamics in a mathematical model of human ventricular tissue with myocytes and fibroblasts.

Cardiac fibroblasts, when coupled functionally with myocytes, can modulate the electrophysiological properties of cardiac tissue. We present systematic numerical studies of such modulation of electrophysiological properties in mathematical models for (a) single myocyte-fibroblast (MF) units and (b) two-dimensional (2D) arrays of such units; our models build on earlier ones and allow for zero-, one-, and two-sided MF couplings. Our studies of MF units elucidate the dependence of the action-potential (AP) morphology on parameters such as [Formula: see text], the fibroblast resting-membrane potential, the fibroblast conductance [Formula: see text], and the MF gap-junctional coupling [Formula: see text]. Furthermore, we find that our MF composite can show autorhythmic and oscillatory behaviors in addition to an excitable response. Our 2D studies use (a) both homogeneous and inhomogeneous distributions of fibroblasts, (b) various ranges for parameters such as [Formula: see text], and [Formula: see text], and (c) intercellular couplings that can be zero-sided, one-sided, and two-sided connections of fibroblasts with myocytes. We show, in particular, that the plane-wave conduction velocity [Formula: see text] decreases as a function of [Formula: see text], for zero-sided and one-sided couplings; however, for two-sided coupling, [Formula: see text] decreases initially and then increases as a function of [Formula: see text], and, eventually, we observe that conduction failure occurs for low values of [Formula: see text]. In our homogeneous studies, we find that the rotation speed and stability of a spiral wave can be controlled either by controlling [Formula: see text] or [Formula: see text]. Our studies with fibroblast inhomogeneities show that a spiral wave can get anchored to a local fibroblast inhomogeneity. We also study the efficacy of a low-amplitude control scheme, which has been suggested for the control of spiral-wave turbulence in mathematical models for cardiac tissue, in our MF model both with and without heterogeneities.


a. A Myocyte-Fibroblast (MF) Composite
The ranges of parameters, which we use for our composite MF system, are consistent with those found in experimental studies and those used in earlier computational studies. For example, in a cell-culture experiment, Rook, et al. [1] have studied rat-heart fibroblasts and reported that the membrane resistance R f , the fibroblast resting membrane potential E f , and the gap-junctional conductance G gap , lie, respectively, in the ranges 3 − 25 GΩ, −20 to −40 mV, and 0.3 − 8.0 nS. Kohl et al. [2] have studied non-excitable cardiac, mechanosensitive fibroblasts from the region of the sinoatrial node in a rat heart. Their study, which uses both intact tissue and cell cultures, estimates that R f ≃ 1 GΩ, E f ≃ −15 ± 10 mV, and G gap ≃ 4 − 6 nS for a well-coupled MF pair. In vitro studies, by Kiseleva, et al. [3], have examined rat mechanosensitive fibroblasts attached to the right atrium; they have found E f ≃ −22 ± 1.9 mV and R f ≃ 0.51 ± 0.01 GΩ, for a control case, and E f ≃ −46.5 ± 1.8 mV and R f ≃ 3.8 ± 0.03 GΩ, in the case of a large infarct caused by a myocardial infarction. In vitro studies by Kamkin, et al. [4] of non-excitable, mechanosensitive, cardiac fibroblasts from the atrium of a human heart have reported E f ≃ −15.9 ± 2.1 mV R f ≃ 4.1 ± 0.1 GΩ. In vitro studies, by Kamkin et al. [5], of rat atrial fibroblasts attached to the sinoatrial node region yielded E f ≃ −22 ± 2 mV and R f ≃ 510 ± 10 M Ω, for the control case, and for the case with myocardial infarction, and E f ≃ −41 ± 3 mV to ≃ −28 ± 3 mV. Recent experiment, in culture, by Chilton et al. [6] have measured the cellular capacitance C f,tot of rat-ventricular fibroblasts by using a patch-clamp recording and found C f,tot ≃ 6.3±1.7 pF; they have shown that the input resistance of fibroblasts R f ≃ 10.7 ± 2.3 GΩ. Their measurements have shown that E f depends on the inwardly rectifying K + current (Kir) and the potassium ion concentration [K + ] o ; e.g, when Kir is expressed, E f is ≃ −65 ± 5 mV and ≃ −80 ± 1.8 mV for [K + ] o = 10 mM and 5.4 mM, respectively. However, when Kir is absent, E f is ≃ −34 ± 2 mV. Furthermore, in culture, Shibukawa, et al. [7] have found, in patch-clamp recordings from rat-ventricular fibroblasts (active), that C f,tot ≃ 4.5 ± 0.4 pF, E f ≃ −58 ± 3.9 mV, R f ≃ 5.5 ± 0.6 GΩ.
The computational studies of mathematical models for fibroblasts, discussed in the "Introduction" Section of the main paper, have also used a wide range of values for parameters for the cellular capacitance C f,tot , the membrane conductance G f , the fibroblast resting membrane potential E f , and the gapjunctional coupling G gap between myocyte and fibroblasts. For example, Xie, et al. [8] have used C f,tot = 25 pF, G f = 0.1 − 4 nS, E f = −50 − 0 mV, and G gap = 0 − 20 nS for an MF composite. The study of Sachse, et al. [9] has used C f,tot = 4.5 pF, E f = −58 mV, and G gap = 0.1 − 100 nS for an MF composite with active fibroblasts. Jacquemet, et al. [10] have studied the MF composite with active fibroblasts by using C f,tot = 4.5 pF, E f = −58 mV, and G gap = 0.09 − 4.05 nS. MacCannell, et al. [11] have used C f,tot = 6 − 60 pF, E f = −49.6 mV and G gap = 1 − 3 nS for their studies of an active-fibroblast model. To investigate in detail the effect of fibroblasts on a myocyte, we use the following wide ranges of parameters (these encompass the ranges used in the experimental and computational studies mentioned above): C f,tot = 6 − 60 pF, G f = 0.1 − 4 nS, E f = −39 to 0 mV, and G gap = 0.3 − 8.0 nS for our MF composites. However, to observe some special properties, such as autorhythmicity of MF composites, we vary the fibroblast parameters and gap-junctional conductances.
It has been noted in Refs. [12][13][14], that a myocyte cell can display autorhythmicity when it is coupled with fibroblasts; in particular, Ref. [12] shows that the cycle length of autorhythmicity activation depends on E f and G gap . We find that G f and C f,tot play a less important role than N f , E f , and G gap in determining whether such autorhythmicity is obtained. In Fig. S4 we give some illustrative plots for N f = 1, E f = 0mV, and G f = 8nS that yield autorhytmicity; Fig. S4 (a) shows a plot of V m versus time t; Fig. S4 (b) contains a plot of the frequency of autorhthymicity f versus G gap for our MF composite; for more detailed studies of the dependence of such autorhythmicity on N f and E f we refer the reader to Ref. [15]. Figure S4 (b) shows that, for the range 0 nS ≲ G gap ≲ 16 nS, the myocyte behaves like an excitable cell, which produces one action potential when it is stimulated electrically; in the range 16 nS ≲ G gap ≲ 23 nS, the myocyte displays autorhthymicity and the cycle length λ cl , the time difference between the upstrokes of two successive action potentials, decreases with increasing G gap ; for G gap ≳ 23nS, the myocyte displays oscillatory behavior. Such autorhythmic and oscillatory responses of an MF composite [15] can occur at lower values of G gap , e.g., G gap = 8 nS, if we increase N f .

b. Wave dynamics in a 2D simulation domain with MF composites
For the case of zero-sided coupling, with E f = 0 mV and G f = 8 nS, Figs. S6(a)-(c) show, respectively, pseudocolor plots of the myocyte transmembrane potential V m , at time t = 2 s, for low-frequency autorhythmicity (e.g., with G gap = 17 nS), high-frequency autorhythmicity (e.g., with G gap = 20 nS), and when the MF composite displays (cf. Fig. S4) oscillatory behavior (e.g., with G gap = 23 nS). In Fig. S6(d), we show the time series of V m (x, y, t), in the interval 0 s ≤ t ≤ 4 s, obtained from three representative points, shown by asterisks in Fig. S6(a), namely, (x = 22.5 mm, y = 22.5 mm) (black filled circles), (x = 67.5 mm, y = 67.5 mm) (blue filled diamonds), and (x = 112.5 mm, 112.5 mm) (red filled triangles); Fig. S6(g) shows the corresponding power spectra, which we calculate from these time series, each of which has 2 × 10 5 data points; each one of these power spectra has discrete, sharp peaks at a fundamental frequency and at its harmonics; the periodic nature of the time series and these peaks in the power spectra provide evidence for the temporally periodic motion of the spiral wave in the region of low-frequency autorhythmicity. The analogs of Figs. S6(d) and (g) are shown in Figs. S6(e) and (h) for the regime of high-frequency autorhythmicity; the time series of V m (x, y, t), in the interval 0 s ≤ t ≤ 4 s, from (x = 22.5 mm, y = 22.5 mm) (black filled circles) and (x = 67.5 mm, y = 67.5 mm) (blue filled diamonds), show irregular behaviors and the corresponding power spectra show subsidiary peaks near the main peaks; however, the time series recorded from (x = 112.5 mm, 112.5 mm) (red filled triangles) shows periodic behavior and, consequently, the corresponding power spectrum has discrete, strong peaks. The analogs of Figs. S6(d) and (g), for the oscillatory regime, are shown in Figs. S6(f) and (i).
If the value of G mm G mf is such that conduction failure occurs in a homogeneous, MF-composite simulation domain, then the MF-composite inhomogeneity behaves somewhat like a conduction inhomogeneity inasmuch as the spiral wave does not enter significantly into the region of the inhomogeneity. To check how far the wave penetrates into the MF-composite inhomogeneity, we show in Figs. S8 (a), (b), and (c) pseudocolor plots of V m at times t = 2s, t = 6s, and t = 8s, respectively, when a square MF-composite inhomogeneity of side ℓ = 33.75mm, with G mm G f f = 1 and G mm G mf = 1 is placed with its bottom left corner at (x = 56.25 mm, y = 56.25 mm). Data for the time series of V m (x, y, t) are recorded at three points of the simulation domain, namely, (x = 90mm, y = 112.5mm), which lies outside the inhomogeneity, (x = 90mm, y = 90mm), at the top-right corner of the inhomogeneity, and (x = 90mm, y = 67.5mm), on the right-middle side of the inhomogeneity; these points are indicated by asterisks in Fig. S8 (c) and the data recorded from them are represented, respectively, by black circles, blue diamonds, and red triangles in Figs. S8 (d)-(f). Figure S8 (d) contains plots of the time series of V m (each one of these time series contain 2 × 10 5 data points). Figure S8 (e) shows the corresponding plots of the inter-beat intervals (ibis) versus the beat number n; and the power spectra E(ω), which follow from the time series of V m , are given in Fig. S8 (f).
From the time series of V m (Fig. S8 (d)), we see small-amplitude oscillations in V m if the time series are obtained from points at the side and corner of the MF-composite inhomogeneity; however, if the point lies outside the inhomogeneity, this time series shows a periodic pattern of action potentials. These times series and the plots of the ibi (Fig. S8 (e)) show that the oscillations in V m , from these three different points, are in phase; to this extent the MF-composite inhomogeneity acts like a conduction inhomogeneity [16]; however, the spiral wave does penetrate the region of the inhomogeneity marginally, so, in this sense, the MF-composite inhomogeneity acts like an ionic inhomogeneity [16].
We turn now to an examination of the interaction of spiral waves with an MF-composite inhomogeneity for different values of G gap . For the same MF-composite inhomogeneity and parameters as in Fig. 13 (c) in the "Results" Section of the main paper, we show in Figs    . Initiation of a spiral wave by the S1-S2 parallel field protocal: To inject a spiral wave, the diffusion constant is set to D mm = 0.000385 cm 2 ms; this is a quarter of its original value, which is 0.00154 cm 2 ms; an S1 stimulus of strength 150 pA/pF is applied for 3 ms at the left boundary; after 560 ms, an S2 stimulus of strength 450 pA/pF is applied for 3 ms just behind the refractory tail of the plane wave initiated by the S1 stimulus; this S2 stimulus is applied over the region x = 360 and 1 ≤ y ≤ 550. We reset D mm to its original value after 880 ms; this procedure yields a fully developed spiral wave at t = 976 ms. The spiral wave configuration in (c) is used as the initial condition for our subsequent studies.    , and (f) for G gap = 0.5 nS, 2 nS, and 8 nS, respectively (data from the points outside and inside the inhomogeneity are represented, respectively, by black circles and red triangles); (g), (h), and (i) show the corresponding plots of the ibi versus the beat number n; and the associated power spectra E(ω) are depicted in (j), (k), and (l).  Table S1. The values of C f,tot , G f , G gap , E f , and N f for a single MF composite and the changes in the AP morphology, relative to that of an uncoupled myocyte. We concentrate on the APD,V max , and V rest and list the changes, indicated by ∆, in these parameters. ∆AP D 70% , ∆AP D 80% , and ∆AP D 90% denote, respectively, the changes in the APD at 70%, 80%, and 90% repolarization. Note that here we have low values (see text) for both G f (0.1 nS) and G gap (0.3 nS).