Ca2+ Cycling Impairment in Heart Failure Is Exacerbated by Fibrosis: Insights Gained From Mechanistic Simulations

Heart failure (HF) is characterized by altered Ca²⁺ cycling, resulting in cardiac contractile dysfunction. Failing myocytes undergo electrophysiological remodeling, which is known to be the main cause of abnormal Ca²⁺ homeostasis. However, structural remodeling, specifically proliferating fibroblasts coupled to myocytes in the failing heart, could also contribute to Ca²⁺ cycling impairment. The goal of the present study was to systematically analyze the mechanisms by which myocyte–fibroblast coupling could affect Ca²⁺ dynamics in normal conditions and in HF. Simulations of healthy and failing human myocytes were performed using established mathematical models, and cells were either isolated or coupled to fibroblasts. Univariate and multivariate sensitivity analyses were performed to quantify effects of ion transport pathways on biomarkers computed from intracellular [Ca²⁺] waveforms. Variability in ion channels and pumps was imposed and populations of models were analyzed to determine effects on Ca²⁺ dynamics. Our results suggest that both univariate and multivariate sensitivity analyses are valuable methodologies to shed light into the ionic mechanisms underlying Ca²⁺ impairment in HF, although differences between the two methodologies are observed at high parameter variability. These can result from either the fact that multivariate analyses take into account ion channels or non-linear effects of ion transport pathways on Ca²⁺ dynamics. Coupling either healthy or failing myocytes to fibroblasts decreased Ca²⁺ transients due to an indirect sink effect on action potential (AP) and thus on Ca²⁺ related currents. Simulations that investigated restoration of normal physiology in failing myocytes showed that Ca²⁺ cycling can be normalized by increasing SERCA and L-type Ca²⁺ current activity while decreasing Na⁺–Ca²⁺ exchange and SR Ca²⁺ leak. Changes required to normalize APs in failing myocytes depended on whether myocytes were coupled to fibroblasts. In conclusion, univariate and multivariate sensitivity analyses are helpful tools to understand how Ca²⁺ cycling is impaired in HF and how this can be exacerbated by coupling of myocytes to fibroblasts. The design of pharmacological actions to restore normal activity should take into account the degree of fibrosis in the failing heart.

Introduction

Heart failure (HF) is a major public health problem worldwide (Savarese and Lund, 2017) and the development of appropriate therapies to manage HF requires an in-depth knowledge of this syndrome. HF is classified as HF with reduced ejection fraction (HFrEF) or HF with preserved ejection fraction (HFpEF) according to left ventricular systolic function and the type of remodeling (Fukuta and Little, 2007). In HFrEF, the heart is unable to pump blood efficiently due to a reduction in cardiac contractility after diverse cardiovascular diseases. The complexity of the excitation–contraction sequence and the multiscale problem can be approached through mathematical models, which significantly help to gain insight into the underlying mechanisms of cardiac dysfunction and guide future research lines (O’Hara et al., 2011; Trayanova and Chang, 2016).

Contractile dysfunction in HF has been associated with an altered Ca²⁺ handling in myocytes, since Ca²⁺ homeostasis is crucial for cell contraction and relaxation (Alpert et al., 2000; Bers, 2000). Failing myocytes present a diminished intracellular Ca²⁺ transient (CaT) with a slow rise time and a reduced rate of [Ca²⁺]_i removal that prolong CaT duration and elevate the diastolic intracellular Ca²⁺ level (Piacentino et al., 2003; Lou et al., 2012). Prolonged action potential duration (APD) and [Na⁺]_i increase are other of the hallmark electrophysiological abnormalities in HF and all of them result from ion channel remodeling in myocytes, i.e., changes in the expression and function of proteins involved in the electrical activity of cells (Gomez et al., 2014b).

A tissue-level hallmark of HF is increased fibrosis and proliferation of cardiac fibroblasts. Several in vitro studies and mathematical modeling studies have documented that electrical coupling between myocytes and fibroblasts will lead to changes in APD and intracellular Ca²⁺ (Zhan et al., 2014; Li et al., 2017). Experimental evidence suggesting the formation of gap junctions between myocytes and fibroblasts in vitro (Gaudesius et al., 2003) has focused researchers attention on the altered electrophysiological properties of myocytes due to this intercellular coupling to explain conduction abnormalities and reentries (Miragoli et al., 2006; Xie et al., 2009a; Majumder et al., 2012). We have already addressed, in a previous work, the consequences on electrical propagation in cardiac tissue under conditions of HF and fibrosis confirming the vulnerability to reentrant activity (Gomez et al., 2014a,b). While electrical changes, with a cellular origin in action potential (AP) properties, have been widely investigated in the heterocellular coupling (Nguyen et al., 2012; Sridhar et al., 2017), changes in Ca²⁺ dynamics have not been explored in depth. It is important, therefore, to understand the role of fibroblasts in the modulation of Ca²⁺ cycling and to progress in the management of HFrEF, improving mechanical contraction.

Therefore, the goal of the present study was to investigate with computational models the effects of fibroblasts on ion transport mechanisms that regulate Ca²⁺ handling in human failing cardiomyocytes. To understand the complex processes taking place in these cells, we made use of sensitivity analyses. Sensitivity calculation has been commonly used for its predictive value in determining electrophysiological properties with parameter variability (Romero et al., 2011; Trenor et al., 2012; Walmsley et al., 2013; Cummins et al., 2014; Mayourian et al., 2016). As univariate and multivariate sensitivity analyses are widely used (Romero et al., 2009; Sobie, 2009), a comparison of both approaches was an initial objective of this work. Inter-subject variability in electrophysiological properties was considered and reproduced by populations of models. Failing populations, with drug-induced alterations in addition to the natural variability, were useful to identify specific combinations of model parameters that could counteract the effects of HF remodeling and fibroblasts. Our results identify the main targets to improve Ca²⁺ dynamics under the pathological conditions explored, improving cardiac contraction recovery.

Materials and Methods

Cellular Models

All simulations were performed at the cellular level. To study the electrophysiological behavior of cardiac myocytes, we used the most complete undiseased human ventricular AP model, developed by O’Hara et al. (2011) (ORd model), which comprises 15 sarcolemmal currents, as shown in Eq. 1, known as fast Na⁺ current (I_Na), late Na⁺ current (I_NaL), transient outward K⁺ current (I_to), L-type Ca²⁺ current (I_CaL), Na⁺ current through the L-type channel (I_CaNa), K⁺ current through the L-type channel (I_CaK), rapid delayed rectifier K⁺ current (I_Kr), slow delayed rectifier K⁺ current (I_Ks), inward rectifier K⁺ current (I_K1), Na⁺/Ca²⁺ exchange current (I_NCX), Na⁺/ K⁺ ATPase current (I_NaK), background currents (I_Nab, I_Cab, I_Kb), and sarcolemmal Ca²⁺ pump current (I_pCa). A detailed Ca²⁺ dynamics is also formulated in the model. Properties such as conductances determining ionic densities and membrane kinetics can be found in the original work (O’Hara et al., 2011). We introduced slight modifications in sodium current formulation, as reported in our previous work (Mora et al., 2017) and leading to ORdmm model, which can also be found in Supplementary Table S1. To reproduce HFrEF phenotype, specific parameters of the model were modified to represent the downregulation or upregulation of cellular proteins experimentally observed in failing cells. This electrophysiological remodeling involved different ion currents and Ca²⁺ fluxes and has already been described in previous studies of our group (Gomez et al., 2014b; Mora et al., 2017). Specifically, the time constant of inactivation of I_NaL (τ_hL), I_NaL conductance, the maximal flux of I_NCX and SR Ca²⁺ leak (J_leak), and the fraction of active binding sites of the Ca⁺² calmodulin-dependent protein kinase II (CaMKa) were upregulated, while conductances of I_to and I_K1, the maximal flux of I_NaK, and Ca²⁺ uptake via SERCA pump (J_SERCA), and SR Ca²⁺-dependence of the steady-state activation of ryanodine receptor release (K_rel,Ca) were downregulated. Further details about values and experimental references can be found in Supplementary Table S2.

$$\begin{array}{c}{\text{I}}_{\text{ion}}\text{=}{\text{I}}_{\text{Na}}{\text{+I}}_{\text{NaL}}{\text{+I}}_{\text{to}}{\text{+I}}_{\text{CaL}}{\text{+I}}_{\text{CaNa}}{\text{+I}}_{\text{CaK}}{\text{+I}}_{\text{Kr}}{\text{+I}}_{\text{Ks}}\text{+}{\text{I}}_{\text{K1}}{\text{+I}}_{\text{NaCa}}{\text{+I}}_{\text{NaK}}{\text{+I}}_{\text{Nab}}{\text{+I}}_{\text{Cab}}{\text{+I}}_{\text{Kb}}{\text{+I}}_{\text{pCa}}\left(1\right)\end{array}$$

The increase in fibroblasts density due to aging and/or HF, and their interaction with myocytes were modeled using established cell–cell coupling equations. Fibroblasts were resistively coupled to one myocyte (Eq. 2), and a current (I_gap) flows from one cell to the other through gap junctions, driven by the potential gradient between the myocyte (V_M) and the fibroblast (V_F) membrane potential and regulated by a junctional conductance (G_gap).

$$\begin{array}{c}{\text{I}}_{\text{gap}}\text{=}{\text{G}}_{\text{gap}}·({\text{V}}_{\text{M}}{\text{-V}}_{\text{F}})\left(2\right)\end{array}$$

The electrical activity of fibroblasts was formulated using an active fibroblast model (MacCannell et al., 2007). The fibroblast AP model incorporates four transmembrane currents: a time- and voltage-dependent delayed-rectifier K⁺ current (I_Kr), an inward-rectifying K⁺ current (I_K1), a Na⁺/K⁺ ATPase current (I_NaK), and a background Na⁺ current (I_b,Na). The membrane potential is governed by the following ordinary differential equation:

$$\begin{array}{c}\frac{{\text{dV}}_{\text{F}}}{\text{dt}}\text{=}\text{-}\frac{\text{1}}{{\text{C}}_{\text{F}}}({\text{I}}_{\text{Kr}}\text{+}{\text{I}}_{\text{K1}}\text{+}{\text{I}}_{\text{NaK}}\text{+}{\text{I}}_{\text{b},\text{Na}}\text{-}{\text{I}}_{\text{gap}})\left(3\right)\end{array}$$

The electrotonic interaction with the myocyte is included in the term I_gap. Regarding fibroblast properties, the membrane capacitance (C_F) was set to 6.3 pF and the resting membrane potential (E_F) had a value of -49.6 mV as in MacCannell et al. (2007). Similarly, the differential equation for the membrane potential of the myocyte, taking into account myocyte–fibroblast coupling, is as follows:

$$\begin{array}{c}\frac{{\text{dV}}_{\text{M}}}{\text{dt}}\text{=-}\frac{\text{1}}{{\text{C}}_{\text{M}}}({\text{I}}_{\text{ion}}{\text{+I}}_{\text{stim}}{\text{+I}}_{\text{gap}})\left(4\right)\end{array}$$

Myocyte dimensions are bigger than fibroblasts, which is considered in membrane capacitance differences (C_M = 153.4 pF vs C_F = 6.3 pF). The value of G_gap was set to 3 nS for a single fibroblast, considered the nominal value in MacCannell et al. (2007) model. This value is within the range of 0.3–8 nS measured in cultured myocyte–fibroblast pairs (Rook et al., 1992). In our simulations, an elevated degree of fibrosis was considered by setting a myocyte–fibroblast ratio of 1:5, obtained by increasing G_gap fivefold (Mayourian et al., 2016; Zimik and Pandit, 2016). In the absence of fibrosis, G_gap is 0, and the myocyte is not coupled to fibroblasts.

Both types of myocytes, with and without electrophysiological HF remodeling, were coupled to the same number of fibroblasts with identical properties to simulate the effect of fibrosis on a single myocyte under different conditions. Four basic models were then considered: ORdmm model, ORdmm model with HF remodeling, ORdmm model with five coupled fibroblasts, and HF ORdmm model coupled to five fibroblasts.

Measurement of Biomarkers

To evaluate the electrophysiological activity of myocytes, and particularly Ca²⁺ dynamics, the electrophysiological indicators chosen were APD from maximal upstroke to 90% of repolarization (APD₉₀), calcium transient (CaT) duration from maximal upstroke to 80% of recovery (CaTD₈₀), CaT rise time from 10 to 90% of upstroke (t_10-90), and systolic and diastolic [Ca²⁺]_i values. All biomarkers were calculated under steady state conditions after application of 1000 stimuli at 1 Hz.

Population of Models

The aforementioned models were used as the baseline to generate four populations of 300 different individuals. Inter-subject electrophysiological variability was represented by variating maximal ion current conductances with a scale factor, assuming that there are variations in the number of ion channels in the cell membrane between individuals (Muszkiewicz et al., 2016). The variation in these scale factors is also a way to simulate the effects of drugs on ion channels, as a simplified action of pharmacological compounds inhibiting or enhancing ion currents. We selected and modified 13 key variables of the model, accounting for maximal ionic conductances and fluxes of I_Na, I_NaL, I_to, I_CaL, I_Kr, I_Ks, I_K1, I_NaK, I_NCX, the SR Ca²⁺ uptake via SERCA pump (J_SERCA), the SR Ca²⁺ release flux via RyR (J_rel), the SR Ca²⁺ leak (J_leak), and I_Nab. These parameters were varied with a scaling factor obtained from a log-normal distribution of standard deviation (σ) equal to 0.3. This led to a 95% (±2σ) of parameters varying between 55 and 182% of its control value, representing inter-individual variability and drug-induced effects. A standard deviation of 0.1 was also considered in different sets of populations.

In order to analyze how failing APs and CaTs could be restored to normal physiological ranges, larger populations of models were generated. Latin Hypercube Sampling (LHS) of parameters in the range of ±60% variation was used for this purpose. In this case, only the basic HF model (with and without coupled fibroblasts) was considered. We set limits for APD₉₀ and systolic and diastolic [Ca²⁺]_i values, in which electrical myocyte activity was considered normal. These limits were obtained from the baseline models as APD_90normal ± (APD_90HF - APD_90normal)/2 and after being compared to experimental studies (Drouin et al., 1998; Li et al., 1998; Pereon et al., 2000). Table 1 shows the considered ranges.

TABLE 1

TABLE 1. Biomarkers ranges for normal electrophysiological activity.

Sensitivity Analysis

Linear regression was employed to analyze biomarker sensitivities to electrophysiological variables. First, a univariate sensitivity analysis was conducted varying each parameter individually by a scaling factor of equal magnitude from the baseline models. The individual variation was set to ±60 or ±15% in different sensitivity analyses. This approach was based on our previous work (Mora et al., 2017).

Multivariate regression analyses were then performed using the generated populations of models. We related all the conductance scaling factors configurations (x) to each biomarker (y) with regression coefficients (b) that attempted to predict the indicators (Eq. 5). These coefficients were obtained by applying partial least squares (PLS) method to data as shown in Eq. 6 and after several transformations, such as log-transformations and z-score calculations (Sobie, 2009; Morotti et al., 2017). The obtained coefficients indicate the relative contribution of parameters variation in biomarker changes and can be considered as sensitivities. This methodology takes into account inter-variable effects that cannot be differentiated in the univariate analyses. Coefficients of determination (R²) indicated the predictive power of the simplified linear model, being R² = 1 a good fit.

$$\begin{array}{c}{\text{y}}_{\text{predicted}}\text{=}{\displaystyle \sum \text{b}}·\text{x}\left(5\right)\end{array}$$

$$\begin{array}{c}{\text{B}}_{\text{PLS}}\text{=}{({\text{X}}^{\text{T}}\text{X})}^{\text{-1}}\times {\text{X}}^{\text{T}}\times \text{Y}\left(6\right)\end{array}$$

Both types of sensitivity analyses were compared with normalized sensitivities to study the similarities and differences in parameter identification as the most important contributors to each electrophysiological characteristic. Normalized sensitivities were calculated as the ratio between each sensitivity and the maximum absolute value obtained for a particular biomarker.

The comparison between isolated and coupled myocytes was performed after adjusting the bar graphs to the standard deviation of the log-transformed biomarkers in N and HF conditions, to ensure equivalent percentages of change.

Results

Comparison of Univariate and Multivariate Sensitivity Analyses in Normal and Failing Conditions

Figure 1 illustrates two populations of models (n = 300) for endocardial human cells paced at 1 Hz, under normal (blue traces) and HF (red traces) conditions, with baseline models for the two respective populations indicated with black solid and dashed lines. The model population calibration process, meant to reproduce natural variability, generates a wide range of physiological AP waveforms and CaTs in both groups. This variability allows us to understand predicted drug effects by considering the HF cells with behavior most similar to normal cells, or vice versa. The electrophysiological remodeling applied to the basic ORdmm model to generate a baseline HF model reproduced the characteristic HF phenotype of prolonged APD and a slower Ca²⁺ dynamics with depressed systolic [Ca²⁺]_i and elevated diastolic levels. In the generated normal and failing populations, significant differences between both conditions are observed in APs and CaTs, although some traces overlap. To understand the variability within the populations and predict what alterations might cause HF cells to behave more like normal cells, we performed a sensitivity analysis.

FIGURE 1

FIGURE 1. Population of normal (N, blue) and heart failure (HF, red) human ventricular cell models. Simulated action potentials (A) and intracellular Ca²⁺ transients (B) obtained from a multivariate set of ionic conductances (n = 300), using the baseline ORdmm model in N conditions (black solid line) and with HF ionic remodeling (black dashed line).

Figure 2 represents the sensitivity of quantitative indicators of AP and CaT characteristics to changes in ion channels obtained in the univariate and multivariate sensitivity analyses. Sensitivities helped reveal the most important modulators of a determined electrophysiological property. High regression coefficients highlighted an important effect of a particular transport mechanism on Ca²⁺ indicators, taking into account the synergy between variables in the case of multivariate analysis. Our results showed that I_Kr was the parameter with the major impact on APD₉₀ (Figure 2A), while SERCA had the leading role in Ca²⁺ biomarkers (Figures 2B–D). Sensitivities were normalized to allow comparison between the sensitivities obtained through a ± 60% univariate sensitivity analysis (light blue for N conditions and pink bars for HF conditions) with multivariate sensitivity analysis (dark blue for N conditions and red bars for HF conditions). In this comparison, results from both methodologies were consistent for APD sensitivities (Figure 2A). Small differences were though observed in the case of Ca²⁺ biomarkers (Figures 2B–D). For instance, Figure 2B shows that in HF, I_Kr seems to have a moderate impact on CaTD₈₀ according to the univariate results (pink bars), while multivariate sensitivity indicates that the effect is lower (red bars). In other cases, with the univariate analysis, the relative importance of parameters can be altered between the N and HF conditions. This was the case for I_CaL, I_NCX, and J_leak on t_10-90 (Figure 2C light blue and pink bars). Figure 2D shows how NCX has the main role in modulating systolic Ca²⁺ according to the univariate sensitivity analysis (light blue in N and pink bars in HF), while the multivariate sensitivity highlights SERCA and I_CaL as the main modulators in N (dark blue bars) and HF (red bar) conditions, respectively, and the exchanger has a secondary role.

FIGURE 2

FIGURE 2. Comparison of univariate and multivariate relative sensitivities in normal (N) and heart failure (HF) conditions. Modulators of (A) action potential duration to 90% of repolarization (APD₉₀), (B) Ca²⁺transient (CaT) duration to 80% of recovery (CaTD₈₀), (C) CaT rise time (t_10-90), and (D) Systolic peak of CaT. Parameter variability of 60%.

Multivariate and univariate analyses also identified similar changes in parameter sensitivities from N to HF (Figure 2, light and dark blue vs pink and red), such as an increased impact of I_NaL on APD₉₀ in HF with respect to normal conditions (Figure 2A), a higher I_Kr influence on CaTD₈₀ (Figure 2B), and a decrease of SERCA modulation effect on systolic [Ca²⁺]_i (Figure 2D).

We further investigated the aforementioned discrepancies by performing additional sensitivity analyses with a lower conductance variability. Parameter variability was reduced to 15% in the univariate study and the standard deviation was decreased to 0.1 in the multivariate regression, resulting in sensitivities slightly different from those described above with higher variability (60% and σ = 0.3, respectively). After reducing parameter variability, negligible differences were found when comparing both methods (see Supplementary Figures S1, S2, light-colored bars). For example, considering low variability in HF (Supplementary Figure S2D, light bars), both analyses (univariate and multivariate) highlighted I_CaL as the most important parameter regulating systolic [Ca²⁺]_i, while univariate sensitivities obtained with 60% variability (darker pink bar) highlighted I_NCX as the main regulator and underestimated the impact of I_CaL.

Coefficients of determination of the multivariable regression are shown in Table 2 for the different sensitivity analyses. Values were closer to 1 as variability decreased. CaTD₈₀ and systolic [Ca²⁺]_i linear fit significantly improved with lower variability, especially in HF (0.938 vs 0.677 and 0.970 vs 0.792, respectively). This indicates that, with large parameter variability, and in the HF condition, significant non-linear relationships between parameter values and Ca²⁺ handling processes make the multivariable regression model less accurate.

TABLE 2

TABLE 2. Coefficients of determination (R²) of the multivariable regression analyses in normal (N) and heart failure (HF).

In summary, univariate and multivariate sensitivity analyses yield very similar results for low variability of parameters, which would correspond to natural inter-subject electrophysiological differences. Univariate analysis is less computational expensive and is thus a valid methodology under these conditions. However, when large variability is applied, which would respond to effects of drugs, significant discrepancies arise between both methodologies. Multivariate analysis should be more reliable, at least for biomarkers with coefficients of determination close to one. Indeed, multivariate analysis considers interaction between the different parameters of the ion channels altered.

Effects of Fibroblast–Myocyte Coupling in Normal and Failing Conditions

When simulations were performed to electrically coupled fibroblasts to a myocyte, the myocyte’s AP and CaT were significantly changed. Figure 3 shows APD and systolic [Ca²⁺]_i reduction in the normal baseline human endocardial model when fibroblasts were coupled (solid vs discontinuous blue trace), as well as in the failing baseline endocardial model (solid vs discontinuous red trace).

FIGURE 3

FIGURE 3. Effect of fibrosis on action potential (A) and Ca²⁺ transient (B). Normal (N) and failing (HF) baseline myocyte models in an uncoupled myocyte (M) and a myocyte coupled to five fibroblasts (M-Fb).

When fibroblasts were coupled to myocytes from two populations of models (n = 300) generated for both N and HF conditions, the above-mentioned effects on APD and CaT were maintained, as shown in Figure 4. CaT traces in normal conditions (Figure 4A) with coupled fibroblasts (discontinuous blue traces) overlapped with CaT traces in the absence of fibroblasts (solid blue traces).

FIGURE 4

FIGURE 4. Action potential and Ca²⁺ transient samples (n = 30) of uncoupled myocytes (solid traces) and myocytes coupled to five fibroblasts (discontinuous traces) from two population of models, obtained from a multivariate set of ionic conductances (n = 300), using the basic ORdmm in normal conditions (A,C) and the basic HF remodeling model (B,D).

To evaluate whether fibroblasts coupling would significantly change sensitivities of electrophysiological biomarkers to parameters variability, multivariate analyses were conducted on these new populations. When a high variability of parameters was considered (σ = 0.3), our results showed that in the presence of fibroblasts, biomarkers sensitivities to ionic variables slightly changed with respect to uncoupled myocytes. As displayed in Figure 5, most of these differences were quantitative, as fibroblasts reduced sensitivities to parameters, but the qualitative role of each parameter was maintained. For instance, SERCA contribution to Ca²⁺ indicators (Figures 5B–D) decreased both in normal and failing conditions when fibroblasts were considered (discontinuous bars).

FIGURE 5

FIGURE 5. Comparison of sensitivities obtained from 4 multivariable regression analyses: normal (N) and heart failure (HF) conditions with or without coupled fibroblasts (Fb). Modulators of (A) action potential duration to 90% of repolarization (APD₉₀), (B) Ca²⁺ transient (CaT) duration to 80% of recovery (CaTD₈₀), (C) CaT rise time (t_10-90), and (D) Systolic peak of CaT. High parameter variability. Regression coefficients (b) are scaled to the standard deviation (σ) of log-normal distributed biomarkers in uncoupled myocytes (M).

To confirm these results, we decreased parameter variability to σ = 0.1 and obtained more accurate regression coefficients. In the new range of parameter variability (see Supplementary Figure S3), sensitivities hardly changed with respect to higher variability.

AP and CaT Restoration in HF

To analyze how failing APs and CaTs could be restored and brought to normal ranges, a larger population of failing models was generated (n = 10,000). The population of models approach generates cells exhibiting a wide variety of physiological behaviors, including cells within the HF population that have APs and CaTs similar to those observed in healthy cells. Selecting this subpopulation and examining the distributions of parameters provides guidance for therapeutic targets. We identified parameter combinations restoring APs and CaTs to waveforms within a healthy range shown in Table 1 in Section “Materials and Methods.” In Figure 6A, red traces represent all models generated from the baseline ORdmm HF model (black dashed line) and blue lines are models satisfying healthy ranges of both AP and CaT biomarkers. Electrophysiological HF phenotype was restored in around 500 models when I_Kr, SERCA, and I_CaL activities were enhanced, and NCX function decreased. Interestingly, NCX distribution presented several outliers at high scaling factors. When myocytes were coupled to fibroblasts (Figure 6B), the role of SERCA, I_CaL, and NCX was similar, whereas I_Kr modulation was not important.

FIGURE 6

FIGURE 6. Restoration of normal action potential (AP) and calcium transient (CaT) in a population of failing models (n = 10,000). All simulated models (red traces) are obtained from a HF basic model (black dashed line) in which parameters have been varied ±60%. Calibrated models (blue traces) are within the limits of normal biomarkers (basic ORdmm in black solid line). Distributions of parameter scaling factors for the restored models. (A) Population of an uncoupled myocyte. (B) Population of a myocyte coupled to five fibroblasts.

Mechanistic Analysis of Calcium Dynamics Impairment in Fibrosis

Fibroblast coupling decreased systolic Ca²⁺ in normal myocytes and exacerbated calcium dynamics impairment in failing myocytes, although parameter sensitivities were hardly affected. To understand the underlying mechanisms, ion currents were carefully analyzed in our simulations. Figure 7 shows the traces of selected electrophysiological variables of the model in the time course of an AP at steady state for a failing myocyte with (solid line) and without (dashed line) fibroblast coupling. To understand how these steady-state conditions were reached, we analyzed transient changes in myocyte [Ca²⁺] after fibroblast coupling, as shown in Figure 8. In these panels, isolated myocytes are initially at steady state, then fibroblasts are coupled beginning on the second beat. As shown in Figure 7C, I_gap is an outward current flowing from the myocyte to fibroblasts, non-existing in the uncoupled myocyte. This outward current decreases membrane potential and has an indirect effect on voltage dependent currents such as I_CaL. Indeed, we can observe that when fibroblasts are coupled (dashed traces) all currents and fluxes present a shorter duration, as well as AP, and a reduced peak (see, for instance, CaT in Figure 7B), except I_CaL. Despite the reduced duration of this current, the initial I_CaL peak contributed to a larger Ca²⁺ influx through these channels (increase from 140 to 146 pC/μF). The integral of NCX was also computed and indicated an increase in inward NCX extruding more Ca²⁺ when fibroblasts were coupled (from 77 to 81 pC/μF). In the transient evolution shown in Figure 8F, the peak of I_CaL exhibited a maintained increase compared with the uncoupled myocyte (shown in solid red). Due to the increased Ca²⁺ extrusion through NCX, however, [Ca²⁺]_JSR progressively declined, leading to the decrease in [Ca²⁺]_i peak observed in steady state.

FIGURE 7

FIGURE 7. Main ionic properties and mechanisms in steady state of a failing myocyte interacting (M-Fb) or not (M) with fibroblasts. (A) Membrane potential (V_m), (B) intracellular Ca²⁺ concentration ([Ca²⁺]_i), (C) sarcoplasmic reticulum (SR) Ca²⁺ concentration ([Ca²⁺]_SR), (D) outward positive current circulating from the myocyte to the fibroblast (I_gap), (E) L-type Ca²⁺ current (I_CaL), (F) Na⁺/Ca²⁺exchange current (I_NCX), (G) SR Ca²⁺ uptake via SERCA pump (J_SERCA), (H) SR Ca²⁺ leak (J_leak), and (I) SR Ca²⁺ release flux via ryanodine receptors (J_rel).

FIGURE 8

FIGURE 8. Transient state of ionic mechanisms and Ca²⁺ concentration in cellular compartments of a failing myocyte coupled to five fibroblasts (dashed lines). Initial steady state of an uncoupled myocyte (solid line). (A) SR Ca²⁺ uptake via SERCA pump (J_SERCA), (B) junctional SR Ca²⁺ concentration ([Ca²⁺]_JSR), (C) SR Ca²⁺ release flux via ryanodine receptors (J_rel), (D) intracellular Ca²⁺ concentration ([Ca²⁺]_i), (E) Na⁺/Ca²⁺exchange current (I_NCX), and (F) L-type Ca²⁺ current (I_CaL).

Discussion

Main Findings

In this study, two methodologies were used to perform sensitivity analyses evaluating the effects of fibroblast–myocyte coupling under normal and HF conditions. Our main findings are that (i) univariate and multivariate analyses yield very similar results and low variability of parameters yields more reliable multivariate sensitivity analyses, (ii) despite exacerbating Ca²⁺ impairment in HF, fibroblast to myocyte coupling does not alter the role of the main mechanisms regulating Ca²⁺ dynamics in myocytes, and (iii) drug action on I_CaL and SERCA enhancement and NCX block would help CaT restoration in HF regardless of fibrosis presence.

Suitable Sensitivity Analyses

Univariate and multivariate sensitivity analyses were performed using an initial large variability (±60%) and then smaller (±15%). At lower variability, changes in electrophysiological properties behaved linearly, making results from univariate and multivariate analysis more similar. Indeed, ionic currents and fluxes work synergistically to generate AP and CaT. Variables such as membrane potential and ionic concentrations link all ionic mechanisms in a way that a change in one parameter (conductivity in this case) has an effect on other parameters, and final differences in electrophysiological properties are the result of an interaction of all these variations. In addition, within the range of variation of a specific parameter, the effects on biomarkers can be different, and assuming linearity can become less accurate when a wide interval is considered. In fact, as electrophysiological remodeling in HF is simulated by changing parameters, the same scale factor applied in N and HF implies a different range of variation of such parameters, accounting for different sensitivities to biomarkers. This explains the increase in APD sensitivity to I_NaL, or the decrease of systolic [Ca²⁺]_i sensitivity to SERCA in HF with respect to N conditions.

Multivariate parameter sets, unlike varying one parameter at a time, can provide information about electrophysiological properties in a wide combination of parameters, which could be useful to evaluate the response of a drug in different individuals, instead of on a fixed baseline model. The high physiological or pathological variability of ionic parameters requires the analysis of the behavior in the whole range.

Despite accuracy loss with high parameter variability, sensitivity analyses are a useful systematic tool to determine the most important mechanisms involved in Ca²⁺ dynamics. We observed that univariate and multivariate sensitivities agreed in the most important parameters contributing to each biomarker. For instance, it was observed the strong impact of I_Kr on APD₉₀, as O’Hara et al. (2011) described in their model. For CaT indicators, the main role of parameters such as SERCA, NCX, and I_CaL was highlighted with both methodologies, agreeing with experimental observations (Ozdemir et al., 2008; Rocchetti et al., 2008; Tamayo et al., 2017), and with a previous study comparing populations of normal and failing myocytes (Walmsley et al., 2013).

In electrophysiological models, biomarker sensitivities are usually calculated within a wide range of variation in ionic conductances, which has been reported experimentally. In their work, Romero et al. (2009) simulated univariate variations up to ±30% to evaluate the arrhythmogenic risk of ionic change, with AP properties falling within the experimental range. In the multivariable regression conducted by Sobie (2009), it was assessed how variability (σ ranging from 0.1 to 0.5) affected the regression model and observed a R² reduction with the increase of σ, although many regression coefficients remained constant, being small the decrease in accuracy obtained for APD. This is consistent with our results, highlighting the robustness of sensitivity analysis for AP biomarkers. However, Ca²⁺-biomarkers sensitivities were more affected by variability. Comparing low and high parameter variability, sensitivities differed more in the univariate case due to changes in parameter–biomarker relation along the variation range, including non-linearities. Sensitivities derived from the multivariate study for HF showed the lowest coefficients of determination, suggesting that HF remodeling enhances non-linearities between variables.

In other modeling studies, experimentally-calibrated populations are first generated with conductances scaled up to ±100%. Constrained biomarkers according to experimental values determine the range of parameters in the subsequent regression analysis (Zhou et al., 2016; Britton et al., 2017). This way, linear fit problem could be solved but natural variability dominates the results and high effects of drugs are not evaluated in this range.

It can be concluded that the roles of the different parameters on electrophysiological biomarkers can be qualitatively estimated using linear methods, including different variability and parameter distributions. Furthermore, our findings highlight that both types of sensitivity studies, univariate and multivariate, provide similar results. However, when large variability is considered, discrepancies arising between these methodologies can become notable, affecting Ca²⁺ properties and patholological conditions to a greater extent. Univariate analysis is less computational expensive and is thus a valid methodology within a moderate variability range and for a reasonable number of parameters. Quantitative results from the multivariate analysis should be more reliable, as parameter interaction is considered, but caution should be taken if coefficients of determination are not close to one, which can happen under HF conditions.

Fibrosis Exacerbates Ca²⁺ Transient Impairment in Heart Failure

This computational study shows that the electrical activity of myocytes, including calcium dynamics, is affected by fibroblast coupling. Specifically, CaT peaks are smaller and SR [Ca²⁺] is reduced when fibroblasts are coupled to myocytes. Despite these alterations, the ionic mechanisms regulating Ca²⁺ cycling are barely affected by the intercellular interaction, indicating that fibroblast coupling does not influence which pathways represent the best drug targets. HF remodeling seems to have a greater impact on the relative role of the different ionic mechanisms that regulate Ca²⁺ cycling. However, our results show that fibroblast coupling could be contributing to the excitation–contraction coupling impairment seen in HF.

An important effect observed when simulating fibroblast–myocyte coupling is the marked APD shortening. Experimental studies showing the existence of gap junctions (connexin43) between fibroblasts and myocytes also revealed changes in AP waveform and conduction velocity due to a modulation of myocyte electrophysiology (Gaudesius et al., 2003; Miragoli et al., 2006; Zlochiver et al., 2008). Previous computational studies exploring fibroblast–myocyte electrophysiological interactions also showed reductions in APD. When MacCannell et al. (2007) developed the mathematical model for the active fibroblast, they coupled it to a human myocyte model and reported an APD shortening, as fibroblasts acted as current sinks. The theoretical work of Jacquemet and Henriquez (2008) also showed that the effect of coupling caused a faster myocyte repolarization, but changing fibroblast properties, such as a less negative fibroblast resting potential, reversed the effect. It is known that fibroblasts differentiate into myofibroblasts in response to inflammation, an activated form which presents contractile proteins, implying the existence of Ca²⁺ cycling (Chilton et al., 2007). To date, no specific models for myofibroblasts have been developed but simulations have been performed increasing the membrane capacitance and depolarizing the resting membrane potential of the fibroblast model (Nguyen et al., 2012; Sridhar et al., 2017). In the present study, we used the fibroblast model by MacCannell et al. (2007) with a resting potential of -49.4 mV, thus the effect in normal and failing myocytes was a reduction in APD. We also conducted some additional simulations (see Supplementary Figure S4) changing to myofibroblast phenotype (C_F = 50 pF, E_F = -24.5 mV), and we observed similar effects, i.e., APD shortening and Ca²⁺-handling impairment. The improvement of myofibroblast models, incorporating Ca²⁺ dynamics, would certainly shed light into the understanding of Ca²⁺ dynamics alterations in the failing tissue.

The number of coupled fibroblasts considered could also alter the results. The uncertainty about the in vivo coupling, different degrees of fibrosis, and the difficult to quantify coupled fibroblasts in tissue has led to the exploration of different values of G_gap or a range of myocyte–fibroblast ratios (Jacquemet and Henriquez, 2008; Sachse et al., 2009; Xie et al., 2009a). We compared the effects of coupling one fibroblast to five (results not shown) and we found that the impact on APD and CaT was lower with one fibroblast and increasing the number of fibroblasts to five did not increase the effects fivefold. This finding suggested that the sensitivity of the myocyte to G_gap is not linear, saturating for higher values, as observed in other analyses of G_gap effects when strong coupling was considered (Jacquemet and Henriquez, 2008). We finally used five fibroblasts in order to represent fibroblast proliferation observed in pathological conditions.

In our simulations, the baseline ORdmm model showed a notable Ca²⁺ impairment when fibroblast interaction was considered in both N and HF conditions, but the population of models revealed that not every individual might have an impaired Ca²⁺ handling in a fibrotic heart. However, in HF, even a minor Ca²⁺-handling alteration should be considered important because it worsens contractility dysfunction. Only a few studies have measured intracellular Ca²⁺ in myocytes interacting with fibroblasts. In co-culture experiments used to investigate the crosstalk between both cell types, it was found that fibroblasts from normal hearts increased CaT amplitude, while fibroblasts and myofibroblasts from pressure-overloaded hearts led to a smaller amplitude associated with a reduction in SR Ca²⁺ content (Cartledge et al., 2015). A more recent in vitro study to explore the influence of heterocellular interactions on cardiomyocyte function showed that only adult fibroblasts had significant consequences on the electrical and mechanical function, by prolonging APD and reducing CaT amplitude (Li et al., 2017). In computational studies, only the work by Zhan et al. (2014) has considered the role of fibroblast proliferation in a mathematical model with Ca²⁺, and reported a longer APD and CaT width, modulating cardiac electromechanical behavior. This discrepancy with our results, a prediction of APD prolongation rather than shortening, might be due to the use of a passive fibroblast model with a resting membrane potential set to -20 mV. Indeed, when modeling fibroblast electrophysiological behavior, the values of several parameters are crucial in the outcome, including the number of fibroblasts coupled to each myocyte, the resting potential, and the capacitance of each fibroblast (Chilton et al., 2005), and whether myofibroblast properties are considered (Sridhar et al., 2017).

Our sensitivity analysis showed minimal differences in the mechanisms determining APD and CaT waveform between myocytes coupled or not to fibroblasts. This could simplify therapies, as presence or absence of fibrotic tissue would not alter the treatment. The slightly reduced sensitivities in the presence of fibrosis indicate that targeted ion transport pathways require a higher variation to produce the same percentage of change to improve the considered electrophysiological property. I_Kr is the main modulator of APD, as O’Hara et al. (2011) reported in their sensitivity analysis, and its enhancement in HF would restore the prolonged APD. However, in myocytes coupled to fibroblasts there is no need to reduce APD via I_Kr modulation as fibroblasts act as a current sink and lead to an earlier repolarization. Regarding CaT indicators, HF remodeling seems to be the main cause of Ca²⁺ alteration, which explains why coupled and uncoupled myocytes require the same measures to restore Ca²⁺ cycling: increase of SERCA and I_CaL, and NCX reduction.

Finally, the mechanistic analysis of Ca²⁺ cycling with fibroblasts reveals that, despite an increase in I_CaL peak with fibroblast coupling, increased NCX Ca²⁺ extrusion leads lower SR content and smaller CaT peaks. In the steady state, there is a balance between fluxes, maintaining constant Ca²⁺ levels, which does not happen in the transient state. In fact, the progressive reduction in intracellular Ca²⁺ is a sign of negative balance between Ca²⁺ influx and efflux in the myocyte, which occurs after a perturbation in the system, i.e., fibroblast coupling. The modulation of all ionic mechanisms including Ca²⁺ cycling is caused by the active role of fibroblasts, becoming a current sink. The new outward current in myocytes, I_gap, accelerated AP repolarization and consequently, all voltage-dependent currents were affected, including I_CaL. The model used for human endocardial AP does not present a notch with a marked early repolarization phase as in other species. For this reason, I_gap does not affect the excitation–contraction coupling by changing the rate of the initial repolarization as in Xie et al. (2009b). In the present study, the implicated mechanism becomes more important in later phases of the AP. We also conducted simulations using the epicardial model of ventricular AP formulated by O’Hara et al. (2011) to see the effect on the early repolarization phase which is indeed present in the epicardial model. As shown in Supplementary Figure S5, no important changes are observed in this phase under the effect of coupled fibroblasts.

The static formulation of myocyte–fibroblast coupling differed from the dynamic gap junctions channels modeled by Brown et al. (2016). They observed a reduction in the junctional current during the upstroke of the AP when considering time- and voltage-sensitive gating channels in homotypic and heterotypic channels, i.e., with different connexins combinations (Cx43 and Cx45). Although it did not significantly alter conduction velocity compared to static gap junctions, fibroblasts with a smaller sink impact could have a different effect on Ca²⁺ dynamics.

Our findings showed that a higher I_CaL peak introduced more Ca²⁺ in the myocyte but did not trigger a higher SR Ca²⁺ release as expected. According to Shannon et al. (2000), there cannot be a release from the SR with a Ca²⁺ load of less than 50% of its maximal content, explaining that the reduced SR Ca²⁺ content characteristic of failing myocytes contributes to reduce force development. Similarly, changes in SR Ca²⁺ load of 58% have been measured in failing isolated myocytes, and related to a smaller CaT (Piacentino et al., 2003). Therefore, a reduction in SR Ca²⁺ content due to fibroblasts could be the cause of an exacerbated Ca²⁺ impairment. The transient evolution shown in Figure 8 helps understand the mechanisms leading to cellular Ca²⁺ loss. One of the advantages of mathematical models is the power to analyze hypothetical situations that cannot be measured experimentally, such as myocyte response to a sudden electrical connection with fibroblasts. In the transient state (Figure 8), we observed that I_CaL peak increased since the first beat, due to the immediate effects of I_gap on AP. Initially, there was also an intracellular Ca²⁺ rise, as SR Ca²⁺ content was still elevated. The mechanisms to remove Ca²⁺ from the cytosol were then activated: SERCA and inward current through NCX (extruding Ca²⁺). NCX role, extruding Ca²⁺ out of the cell, becomes relevant in the transient state because if there is an imbalance with Ca²⁺ influx through I_CaL channels, Ca²⁺ loss occurs. This mechanism might explain the reduced SR Ca²⁺ load in myocytes interacting with fibroblasts.

Our findings suggest that fibroblasts increase Ca²⁺ impairment in HF by further reducing SR Ca²⁺ content in myocytes.

Limitations

Several limitations need to be considered when drawing conclusions from the present study. Although simulations of cellular electrophysiological behavior and systematic analyses of biomarkers complement and enrich experimental research, uncertainties in the development of the mathematical models employed might affect the outcome of the simulations. The limited availability of electrophysiology data from human fibroblasts has led to the use of a fibroblast model developed on the basis of adult rat ventricular tissue data, without taking into consideration potential changes in HF. Fibroblasts with the same characteristics were considered in normal and failing conditions, but electrophysiological remodeling in fibroblast currents could have an additional effect on myocytes (Aguilar et al., 2014). Another limitation is the use of a fibroblast model, instead of using a myofibroblast model (not available to date), which is the characteristic form in pathological conditions, which would take into account Ca²⁺ dynamics in these cells. Ca²⁺ signaling has been observed in human cardiac fibroblasts although with different pathways to those in contractile myocytes (Chen et al., 2010), and myofibroblasts show intracellular CaTs, modulated by intercellular coupling with myocytes (Chilton et al., 2007), but neither of them have been modeled. Specific data for fibroblasts features are still lacking and one of the major concerns in studies focusing myocyte–fibroblast interactions is the use of data from co-cultures and the need of data describing the real electrophysiological behavior in vivo.

Gap junctional coupling between a myocyte and a fibroblast has been modeled by a simple conductance as in previous simulation studies of fibroblast-myocyte coupling. However, recent findings suggesting the interaction of connexin45 with CaM (Zou et al., 2017) highlight that some ions such as Ca²⁺ can modulate the coupling current. If the modulation of the intercellular coupling by Ca²⁺ ions was taken into account, larger alterations in Ca²⁺ dynamics mechanisms could be observed. Any mechanism related to Ca²⁺ included into our model, such as realistic Ca²⁺ dynamics in fibroblasts, may affect the regulation of Ca²⁺ cycling in myocytes, and enhance the influence of fibroblasts.

Other structural modifications in HF, concerning myocytes, involve the loss of transverse tubules and reorganization of the cell membrane (Lyon et al., 2009). We have assumed that with the electrophysiological remodeling applied to ionic mechanisms, we qualitatively reproduce the delayed SR Ca²⁺ release resulting from detubulation. However, a detailed model including local changes in membrane structure related to the spatial organization of ion channels in HF could provide more accurate results about Ca²⁺ impairment (Nivala et al., 2015; Sanchez-Alonso et al., 2016).

Although Ca²⁺ homeostasis is related to contraction and relaxation force, other mechanical factors, such as myofibroblast contractility when fibroblasts are differentiated into the activated form or extracellular stiffness due to the excessive accumulation of collagen in HF, can contribute to cardiac dysfunction. Our model only considers intracellular Ca²⁺ in myocytes at a cellular level to evaluate the contraction of the heart. The behavior in tissue or even in the whole organ could also be different. A defined architecture of the myocyte when considering spatial distribution would allow to take into consideration the distribution of gap junctions as well as the extracellular space which could modulate Ca²⁺ homeostasis as it does in electrical propagation (Cabo and Boyden, 2009; Sachse et al., 2009; Seidel et al., 2010; Greisas and Zlochiver, 2016). We consider a natural continuation of the present work the analysis of fibrosis effects on calcium waves in 2D and 3D tissues, in which altered impulse propagation and generation of arrhythmias should develop because of the heterogeneities in tissue, according to other studies (Zlochiver et al., 2008; Sachse et al., 2009; Nguyen et al., 2012; Gomez et al., 2014a; Greisas and Zlochiver, 2016; Zimik and Pandit, 2016).

Regarding sensitivity analyses, the non-linearities are the main issue in quantifying the impact of parameters on biomarkers. To restore the electrical activity and contraction of myocytes, the effect of some drugs can imply a high change in an ionic mechanism, and the predictive value of this methodology decreases with increased variability.

Nonetheless, despite these limitations, cellular simulations can shed light in the causes of Ca²⁺ impairment observed in HF, which together with the arrhythmogenic activity of fibrotic tissue, can compromise the function of the myocardium. In the present study, the measures suggested to restore Ca²⁺ dynamics and contractility at cellular level are reliable because they are supported by experiments and according to our findings, tackling the electrophysiology remodeling in failing myocytes can also improve the effects of fibroblasts interactions.

Author Contributions

MM, JF, and BT conceived and designed the study, and analyzed the data. MM performed simulations. JG and ES contributed to the methodology. MM and BT wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was partially supported by the National Science Foundation (MCB 1615677), the American Heart Association (15GRNT25490006), the “Plan Estatal de Investigación Científica y Técnica y de Innovación 2013–2016 from the Ministerio de Economía, Industria y Competitividad of Spain and Fondo Europeo de Desarrollo Regional (FEDER) DPI2016-75799-R (AEI/FEDER, UE)”, and the “Programa de Ayudas de Investigación y Desarrollo (PAID-01-17)” from the Universitat Politècnica de València. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.01194/full#supplementary-material

References

Aguilar, M., Qi, X. Y., Huang, H., Comtois, P., and Nattel, S. (2014). Fibroblast electrical remodeling in heart failure and potential effects on atrial fibrillation. Biophys. J. 107, 2444–2455. doi: 10.1016/j.bpj.2014.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Alpert, N. R., Hasenfuss, G., Leavitt, B. J., Ittleman, F. P., Pieske, B., and Mulieri, L. A. (2000). A mechanistic analysis of reduced mechanical performance in human heart failure. Jpn Heart J. 41, 103–115. doi: 10.1536/jhj.41.103

PubMed Abstract | CrossRef Full Text | Google Scholar

Bers, D. M. (2000). Calcium fluxes involved in control of cardiac myocyte contraction. Circ. Res. 87, 275–281. doi: 10.1161/01.RES.87.4.275

CrossRef Full Text | Google Scholar

Britton, O. J., Bueno-Orovio, A., Virág, L., Varró, A., and Rodriguez, B. (2017). The electrogenic Na⁺/K⁺ pump is a key determinant of repolarization abnormality susceptibility in human ventricular cardiomyocytes: a population-based simulation study. Front. Physiol. 8:278. doi: 10.3389/fphys.2017.00278

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, T. R., Krogh-Madsen, T., and Christini, D. J. (2016). Illuminating myocyte-fibroblast homotypic and heterotypic gap junction dynamics using dynamic clamp. Biophys. J. 111, 785–797. doi: 10.1016/j.bpj.2016.06.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Cabo, C., and Boyden, P. A. (2009). Extracellular space attenuates the effect of gap junctional remodeling on wave propagation: a computational study. Biophys. J. 96, 3092–3101. doi: 10.1016/j.bpj.2009.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Cartledge, J. E., Kane, C., Dias, P., Tesfom, M., Clarke, L., Mckee, B., et al. (2015). Functional crosstalk between cardiac fibroblasts and adult cardiomyocytes by soluble mediators. Cardiovasc. Res. 105, 260–270. doi: 10.1093/cvr/cvu264

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J.-B., Tao, R., Sun, H.-Y., Tse, H.-F., Lau, C.-P., and Li, G.-R. (2010). Multiple Ca²⁺ signaling pathways regulate intracellular Ca²⁺ activity in human cardiac fibroblasts. J. Cell. Physiol. 223, 68–75. doi: 10.1002/jcp.22010

PubMed Abstract | CrossRef Full Text | Google Scholar

Chilton, L., Giles, W. R., and Smith, G. L. (2007). Evidence of intercellular coupling between co-cultured adult rabbit ventricular myocytes and myofibroblasts. J. Physiol. 583, 225–236. doi: 10.1113/jphysiol.2007.135038

PubMed Abstract | CrossRef Full Text | Google Scholar

Chilton, L., Ohya, S., Freed, D., George, E., Drobic, V., Shibukawa, Y., et al. (2005). K+ currents regulate the resting membrane potential, proliferation, and contractile responses in ventricular fibroblasts and myofibroblasts. Am. J. Physiol. Circ. Physiol. 288, H2931–H2939. doi: 10.1152/ajpheart.01220.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Cummins, M. A., Dalal, P. J., Bugana, M., Severi, S., and Sobie, E. A. (2014). Comprehensive analyses of ventricular myocyte models identify targets exhibiting favorable rate dependence. PLoS Comput. Biol. 10:e1003543. doi: 10.1371/journal.pcbi.1003543

PubMed Abstract | CrossRef Full Text | Google Scholar

Drouin, E., Lande, G., and Charpentier, F. (1998). Amiodarone reduces transmural heterogeneity of repolarization in the human heart. J. Am. Coll. Cardiol. 32, 1063–1067. doi: 10.1016/S0735-1097(98)00330-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukuta, H., and Little, W. C. (2007). Contribution of systolic and diastolic abnormalities to heart failure with a normal and a reduced ejection fraction. Prog. Cardiovasc. Dis. 49, 229–240. doi: 10.1016/j.pcad.2006.08.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaudesius, G., Miragoli, M., Thomas, S. P., and Rohr, S. (2003). Coupling of cardiac electrical activity over extended distances by fibroblasts of cardiac origin. Circ. Res. 93, 421–428. doi: 10.1161/01.RES.0000089258.40661.0C

PubMed Abstract | CrossRef Full Text | Google Scholar

Gomez, J. F., Cardona, K., Martinez, L., Saiz, J., and Trenor, B. (2014a). Electrophysiological and structural remodeling in heart failure modulate arrhythmogenesis. 2D simulation study. PLoS One 9:e103273. doi: 10.1371/journal.pone.0103273

PubMed Abstract | CrossRef Full Text | Google Scholar

Gomez, J. F., Cardona, K., Romero, L., Ferrero, J. M., and Trenor, B. (2014b). Electrophysiological and structural remodeling in heart failure modulate arrhythmogenesis. 1D simulation study. PLoS One 9:e106602. doi: 10.1371/journal.pone.0106602

PubMed Abstract | CrossRef Full Text | Google Scholar

Greisas, A., and Zlochiver, S. (2016). The multi-domain fibroblast/myocyte coupling in the cardiac tissue: a theoretical study. Cardiovasc. Eng. Technol. 7, 290–304. doi: 10.1007/s13239-016-0266-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacquemet, V., and Henriquez, C. S. (2008). Loading effect of fibroblast-myocyte coupling on resting potential, impulse propagation, and repolarization: insights from a microstructure model. Am. J. Physiol. Circ. Physiol. 294, H2040–H2052. doi: 10.1152/ajpheart.01298.2007

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, G. R., Feng, J., Yue, L., and Carrier, M. (1998). Transmural heterogeneity of action potentials and Ito1 in myocytes isolated from the human right ventricle. Am. J. Physiol. 275, H369–H377.

PubMed Abstract | Google Scholar

Li, Y., Asfour, H., and Bursac, N. (2017). Age-dependent functional crosstalk between cardiac fibroblasts and cardiomyocytes in a 3D engineered cardiac tissue. Acta Biomater. 55, 120–130. doi: 10.1016/j.actbio.2017.04.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Lou, Q., Janks, D. L., Holzem, K. M., Lang, D., Onal, B., Ambrosi, C. M., et al. (2012). Right ventricular arrhythmogenesis in failing human heart: the role of conduction and repolarization remodeling. Am. J. Physiol. Heart Circ. Physiol. 303, H1426–H1434. doi: 10.1152/ajpheart.00457.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyon, A. R., MacLeod, K. T., Zhang, Y., Garcia, E., Kanda, G. K., Lab, M. J., et al. (2009). Loss of T-tubules and other changes to surface topography in ventricular myocytes from failing human and rat heart. Proc. Natl. Acad. Sci. U.S.A. 106, 6854–6859. doi: 10.1073/pnas.0809777106

PubMed Abstract | CrossRef Full Text | Google Scholar

MacCannell, K. A., Bazzazi, H., Chilton, L., Shibukawa, Y., Clark, R. B., and Giles, W. R. (2007). A mathematical model of electrotonic interactions between ventricular myocytes and fibroblasts. Biophys. J. 92, 4121–4132. doi: 10.1529/biophysj.106.101410

CrossRef Full Text | Google Scholar

Majumder, R., Nayak, A. R., and Pandit, R. (2012). Nonequilibrium arrhythmic states and transitions in a mathematical model for diffuse fibrosis in human cardiac tissue. PLoS One 7:e45040. doi: 10.1371/journal.pone.0045040

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayourian, J., Savizky, R. M., Sobie, E. A., and Costa, K. D. (2016). Modeling electrophysiological coupling and fusion between human mesenchymal stem cells and cardiomyocytes. PLoS Comput. Biol. 12:e1005014. doi: 10.1371/journal.pcbi.1005014

PubMed Abstract | CrossRef Full Text | Google Scholar

Miragoli, M., Gaudesius, G., and Rohr, S. (2006). Electrotonic modulation of cardiac impulse conduction by myofibroblasts. Circ. Res. 98, 801–810. doi: 10.1161/01.RES.0000214537.44195.a3

PubMed Abstract | CrossRef Full Text | Google Scholar

Mora, M. T., Ferrero, J. M., Romero, L., and Trenor, B. (2017). Sensitivity analysis revealing the effect of modulating ionic mechanisms on calcium dynamics in simulated human heart failure. PLoS One 12:e0187739. doi: 10.1371/journal.pone.0187739

PubMed Abstract | CrossRef Full Text | Google Scholar

Morotti, S., Nieves-Cintrón, M., Nystoriak, M. A., Navedo, M. F., and Grandi, E. (2017). Predominant contribution of L-type Cav1.2 channel stimulation to impaired intracellular calcium and cerebral artery vasoconstriction in diabetic hyperglycemia. Channels 11, 340–346. doi: 10.1080/19336950.2017.1293220

PubMed Abstract | CrossRef Full Text | Google Scholar

Muszkiewicz, A., Britton, O. J., Gemmell, P., Passini, E., Sánchez, C., Zhou, X., et al. (2016). Variability in cardiac electrophysiology: using experimentally-calibrated populations of models to move beyond the single virtual physiological human paradigm. Prog. Biophys. Mol. Biol. 120, 115–127. doi: 10.1016/j.pbiomolbio.2015.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, T. P., Xie, Y., Garfinkel, A., Qu, Z., and Weiss, J. N. (2012). Arrhythmogenic consequences of myofibroblast–myocyte coupling. Cardiovasc. Res. 93, 242–251. doi: 10.1093/cvr/cvr292

PubMed Abstract | CrossRef Full Text | Google Scholar

Nivala, M., Song, Z., Weiss, J. N., and Qu, Z. (2015). T-tubule disruption promotes calcium alternans in failing ventricular myocytes: mechanistic insights from computational modeling. J. Mol. Cell. Cardiol. 79, 32–41. doi: 10.1016/j.yjmcc.2014.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Hara, T., Virág, L., Varró, A., and Rudy, Y. (2011). Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS Comput. Biol. 7:e1002061. doi: 10.1371/journal.pcbi.1002061

PubMed Abstract | CrossRef Full Text | Google Scholar

Ozdemir, S., Bito, V., Holemans, P., Vinet, L., Mercadier, J. J., Varro, A., et al. (2008). Pharmacological inhibition of Na/Ca exchange results in increased cellular Ca²⁺ load attributable to the predominance of forward mode block. Circ. Res. 102, 1398–1405. doi: 10.1161/CIRCRESAHA.108.173922

PubMed Abstract | CrossRef Full Text | Google Scholar

Pereon, Y., Demolombe, S., Baro, I., Drouin, E., Charpentier, F., and Escande, D. (2000). Differential expression of KvLQT1 isoforms across the human ventricular wall. Am. J. Physiol. Hear. Circ. Physiol. 278, H1908–H1915. doi: 10.1152/ajpheart.2000.278.6.H1908

PubMed Abstract | CrossRef Full Text | Google Scholar

Piacentino, V., Weber, C. R., Chen, X., Weisser-Thomas, J., Margulies, K. B., Bers, D. M., et al. (2003). Cellular basis of abnormal calcium transients of failing human ventricular myocytes. Circ. Res. 92, 651–658. doi: 10.1161/01.RES.0000062469.83985.9B

PubMed Abstract | CrossRef Full Text | Google Scholar

Rocchetti, M., Alemanni, M., Mostacciuolo, G., Barassi, P., Altomare, C., Chisci, R., et al. (2008). Modulation of sarcoplasmic reticulum function by PST2744 [istaroxime; (E,Z)-3-((2-aminoethoxy)imino) androstane-6,17-dione hydrochloride)] in a pressure-overload heart failure model. J. Pharmacol. Exp. Ther. 326, 957–965. doi: 10.1124/jpet.108.138701

PubMed Abstract | CrossRef Full Text | Google Scholar

Romero, L., Carbonell, B., Trenor, B., Rodríguez, B., Saiz, J., and Ferrero, J. M. (2011). Systematic characterization of the ionic basis of rabbit cellular electrophysiology using two ventricular models. Prog. Biophys. Mol. Biol. 107, 60–73. doi: 10.1016/J.PBIOMOLBIO.2011.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Romero, L., Pueyo, E., Fink, M., and Rodríguez, B. (2009). Impact of ionic current variability on human ventricular cellular electrophysiology. Am. J. Physiol. Heart Circ. Physiol. 297, H1436–H1445. doi: 10.1152/ajpheart.00263.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Rook, M. B., van Ginneken, A. C., de Jonge, B., el Aoumari, A., Gros, D., and Jongsma, H. J. (1992). Differences in gap junction channels between cardiac myocytes, fibroblasts, and heterologous pairs. Am. J. Physiol. Physiol. 263, C959–C977. doi: 10.1152/ajpcell.1992.263.5.C959

PubMed Abstract | CrossRef Full Text | Google Scholar

Sachse, F. B., Moreno, A. P., Seemann, G., and Abildskov, J. A. (2009). A model of electrical conduction in cardiac tissue including fibroblasts. Ann. Biomed. Eng. 37, 874–889. doi: 10.1007/s10439-009-9667-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez-Alonso, J. L., Bhargava, A., O’Hara, T., Glukhov, A. V., Schobesberger, S., Bhogal, N., et al. (2016). Microdomain-specific modulation of L-type calcium channels leads to triggered ventricular arrhythmia in heart failure. Circ. Res. 119, 944–945. doi: 10.1161/CIRCRESAHA.116.308698

PubMed Abstract | CrossRef Full Text | Google Scholar

Savarese, G., and Lund, L. H. (2017). Global public health burden of heart failure. Card. Fail. Rev. 3, 7–11. doi: 10.15420/cfr.2016:25:2

CrossRef Full Text | Google Scholar

Seidel, T., Salameh, A., and Dhein, S. (2010). A simulation study of cellular hypertrophy and connexin lateralization in cardiac tissue. Biophys. J. 99, 2821–2830. doi: 10.1016/j.bpj.2010.09.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, T. R., Ginsburg, K. S., and Bers, D. M. (2000). Potentiation of fractional sarcoplasmic reticulum calcium release by total and free intra-sarcoplasmic reticulum calcium concentration. Biophys. J. 78, 334–343. doi: 10.1016/S0006-3495(00)76596-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Sobie, E. A. (2009). Parameter sensitivity analysis in electrophysiological models using multivariable regression. Biophys. J. 96, 1264–1274. doi: 10.1016/j.bpj.2008.10.056

PubMed Abstract | CrossRef Full Text | Google Scholar

Sridhar, S., Vandersickel, N., and Panfilov, A. V. (2017). Effect of myocyte-fibroblast coupling on the onset of pathological dynamics in a model of ventricular tissue. Sci. Rep. 7:40985. doi: 10.1038/srep40985

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamayo, M., Manzanares, E., Bas, M., Martín-Nunes, L., Val-Blasco, A., Jesús Larriba, M., et al. (2017). Calcitriol (1,25-dihydroxyvitamin D3) increases L-type calcium current via protein kinase A signaling and modulates calcium cycling and contractility in isolated mouse ventricular myocytes. Heart Rhythm 14, 432–439. doi: 10.1016/j.hrthm.2016.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Trayanova, N. A., and Chang, K. C. (2016). How computer simulations of the human heart can improve anti-arrhythmia therapy. J. Physiol. 594, 2483–2502. doi: 10.1113/JP270532

PubMed Abstract | CrossRef Full Text | Google Scholar

Trenor, B., Cardona, K., Gomez, J. F., Rajamani, S., Ferrero, J. M., Belardinelli, L., et al. (2012). Simulation and mechanistic investigation of the arrhythmogenic role of the late sodium current in human heart failure. PLoS One 7:e32659. doi: 10.1371/journal.pone.0032659

PubMed Abstract | CrossRef Full Text | Google Scholar

Walmsley, J., Rodriguez, J. F., Mirams, G. R., Burrage, K., Efimov, I. R., and Rodriguez, B. (2013). mRNA expression levels in failing human hearts predict cellular electrophysiological remodeling: a population-based simulation study. PLoS One 8:e56359. doi: 10.1371/journal.pone.0056359

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Garfinkel, A., Camelliti, P., Kohl, P., Weiss, J. N., and Qu, Z. (2009a). Effects of fibroblast-myocyte coupling on cardiac conduction and vulnerability to reentry: a computational study. Heart Rhythm 6, 1641–1649. doi: 10.1016/j.hrthm.2009.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Garfinkel, A., Weiss, J. N., and Qu, Z. (2009b). Cardiac alternans induced by fibroblast-myocyte coupling: mechanistic insights from computational models. Am. J. Physiol. Circ. Physiol. 297, H775–H784. doi: 10.1152/ajpheart.00341.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhan, H., Xia, L., Shou, G., Zang, Y., Liu, F., and Crozier, S. (2014). Fibroblast proliferation alters cardiac excitation conduction and contraction: a computational study. J. Zhejiang Univ. Sci. B 15, 225–242. doi: 10.1631/jzus.B1300156

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Bueno-Orovio, A., Orini, M., Hanson, B., Hayward, M., Taggart, P., et al. (2016). In vivo and in silico investigation into mechanisms of frequency dependence of repolarization alternans in human ventricular cardiomyocytes. Circ. Res. 118, 266–278. doi: 10.1161/CIRCRESAHA.115.307836

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimik, S., and Pandit, R. (2016). Instability of spiral and scroll waves in the presence of a gradient in the fibroblast density: the effects of fibroblast-myocyte coupling. New J. Phys. 18, 1–13. doi: 10.1088/1367-2630/18/12/123014

CrossRef Full Text | Google Scholar

Zlochiver, S., Muñoz, V., Vikstrom, K. L., Taffet, S. M., Berenfeld, O., and Jalife, J. (2008). Electrotonic myofibroblast-to-myocyte coupling increases propensity to reentrant arrhythmias in two-dimensional cardiac monolayers. Biophys. J. 95, 4469–4480. doi: 10.1529/biophysj.108.136473

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, J., Salarian, M., Chen, Y., Zhuo, Y., Brown, N. E., Hepler, J. R., et al. (2017). Direct visualization of interaction between calmodulin and connexin45. Biochem. J. 474, 4035–4051. doi: 10.1042/BCJ20170426

PubMed Abstract | CrossRef Full Text | Google Scholar