We first used a mouse ventricular myocyte model to investigate the effects of changing a single electrophysiological parameter on the CaT. This model (see Materials and methods) combines elements of previously published ventricular mycoyte models (Shiferaw et al., 2003; Shannon et al., 2004; Bondarenko et al., 2004; Mahajan et al., 2008). The CaT was characterized by its amplitude, defined as the difference ${\displaystyle \mathrm{\Delta }{\left[\text{Ca}\right]}_{\text{i}}}$ between the diastolic and peak value of the cytosolic Ca²⁺ concentration ${\left[\text{Ca}\right]}_{\text{i}}$, and the time-averaged value of ${\left[\text{Ca}\right]}_{\text{i}}$ over one pacing period, denoted by ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$. The CaT amplitude ${\displaystyle \mathrm{\Delta }{\left[\text{Ca}\right]}_{\text{i}}}$ is a major determinant of the contractile force while ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$ provides an average measure of the cytosolic Ca²⁺ concentration in the cell. Both quantities will be used as Ca²⁺ sensors for our multi-parameter search of GES and examining individual parameter effects will be useful later to interpret the results of that search. We vary the conductances of sarcolemmal currents and transporters depicted in Figure 2A and the expression levels of Ca²⁺ handling proteins that include the ryanodine receptor (RyR) Ca²⁺ release channels and the sarcoplasmic reticulum (SR) Ca²⁺ ATPase SERCA, which pumps Ca²⁺ from the cytosol into the SR. For each parameter value, we pace the myocyte at a 4 Hz frequency for many beats until a steady-state is reached where the CaT profile used to calculate ${\displaystyle \mathrm{\Delta }{\left[\text{Ca}\right]}_{\text{i}}}$ and ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$ and the intracellular sodium concentration ${\left[\text{Na}\right]}_{\text{i}}$ no longer vary from beat to beat.

Figure 2

Download asset Open asset

Figure 2 shows the effects of individual parameter changes on the steady-state CaT amplitude (Figure 2A) and average ${\left[\text{Ca}\right]}_{\text{i}}$ (Figure 2B). Those six parameters were selected because they control the major currents influencing the CaT. Both quantities are plotted as a function of conductance fold change G/G_ref where G_ref is a reference value producing a normal CaT. Increasing the conductance of the inward L-type Ca²⁺ current ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ is seen to strongly increase both ${\displaystyle \mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$ and ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$ but has a weak effect on the AP waveform (Figure 2C). The effect on the CaT stem from the fact that ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ is the main trigger of Ca²⁺-induced Ca²⁺ release(CICR), which transfers a large amount of Ca²⁺ from the SR to the cytosol. The weak effect on APD is due to the fact that the increase of CaT amplitude causes ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ to inactivate more rapidly, thereby opposing AP prolongation. In contrast, increasing the conductances of K⁺ currents that are dominant in the mouse such as ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$ (fast inactivating component of the transient outward current) and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ causes both ${\displaystyle \mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$ and ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$ to decrease. Increasing either of those K⁺ currents speeds up repolarization (as illustrated for ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$ in Figure 2D) and hence inactivation of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, thereby reducing the magnitude of CICR. Increasing the conductance of the sodium-calcium exchanger current ${\displaystyle {\mathrm{I}}_{\text{NaCa}}}$ also causes both ${\displaystyle \mathrm{\Delta }[\mathrm{C}\mathrm{a}{]}_{\mathrm{i}}}$ and ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$ to decrease by enhancing the forward mode of this current that extrudes Ca²⁺ from the cytosol. The results of Figure 2A are consistent with a study of the effects of single conductance change in rat ventricular myocytes (Devenyi and Sobie, 2016), albeit with a much stronger influence of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ conductance on CaT amplitude in the present mouse model. Changing RyR expression from its reference value is seen to leave ${\displaystyle \mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$ and ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$ almost unchanged, even though it strongly affects the SR Ca²⁺ concentration ${\left[\text{Ca}\right]}_{\text{SR}}$ (Figure 2E). This behavior reflects the well-known effect that making RyR channels more leaky (e.g. by addition of caffeine that increases RyR activity or, similarly here, by increasing the magnitude of the Ca²⁺ release flux through RyRs) yields a transient increase in CaT amplitude, but no change in the steady-state CaT amplitude after ${\left[\text{Ca}\right]}_{\text{SR}}$ adjusts to a lower steady-state level (Bers, 2001). This effect is illustrated by time traces of ${\left[\text{Ca}\right]}_{\text{SR}}$ and ${\left[\text{Ca}\right]}_{\text{i}}$ in steady-state for different RyR expression levels in Figure 2E. Finally, changing the expression level of SERCA has opposite effects on ${\displaystyle \mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$ and ${\displaystyle ⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$. Increasing SERCA magnitude increases SR Ca²⁺ load, thereby increasing the amount of SR Ca²⁺ release and CaT amplitude, but at the same time depletes Ca²⁺ from the cytosol.

Next, we performed a computational search for combinations of parameters that yield a normal electrophysiological output as defined by the steady-state CaT amplitude ${\displaystyle \mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$, time averaged cytosolic Ca²⁺ concentration ${\displaystyle ⟨{\left[Ca\right]}_{\mathrm{i}}⟩}$, and intracellular sodium concentration ${\left[\text{Na}\right]}_{\text{i}}$ at a 4 Hz pacing frequency. A GES search that uses the diastolic and peak ${\left[\text{Ca}\right]}_{\text{i}}$ values as Ca²⁺ sensors, instead of ${\displaystyle \mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$ (the difference between the peak and diastolic ${\left[\text{Ca}\right]}_{\text{i}}$ values) and time averaged ${\left[\text{Ca}\right]}_{\text{i}}$, gives nearly identical results. So our Ca²⁺ sensors can be straightforwardly interpreted physiologically as requirements of normal diastolic and systolic contractile function necessary for a normal arterial blood pressure at the organism scale. A ‘good enough solution’ (GES) was defined as a combination of electrophysiological parameters that produces output values of those three quantities that are close enough to normal target values, which are defined as the values ${\displaystyle \mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}^{\ast }}$, ${\displaystyle ⟨{\left[Ca\right]}_{\mathrm{i}}{⟩}^{\ast }}$, and ${\left[\text{Na}\right]}_{\text{i}}^{\ast }$ corresponding to the reference set of parameters (G_ref values) of the ventricular mycoyte model. The search was conducted by defining a cost function

which is an aggregate measure of the deviation of output sensors $S_n\left(p\right)$ from their desired target values $S_n^{\ast }$. Here, $N=3$ with ${\displaystyle S_1=\mathrm{\Delta }{\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$, ${\displaystyle S_2=⟨{\left[\text{Ca}\right]}_{\text{i}}⟩}$, and $S_3={\left[\text{Na}\right]}_{\text{i}}$, and $ϵ$ is a small tolerance that we choose to be 5%. $E$ is a function of model parameters ${\displaystyle \mathbf{p}=\left(p_1,p_2,\dots \right)}$ chosen to consist of the conductances of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$, ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$, and ${\displaystyle {\mathrm{I}}_{\text{NaCa}}}$ as well as RyR and SERCA expression levels. Effects of individual changes of those parameters on CaT properties measured by $S_1$ and $S_2$ are shown in Figure 2A,B. Conductances of other sarcolemmal currents that were found to have a negligible effect on the CaT were kept constant. The search for GES was conducted by first generating a large population of $\sim 10,000$ randomly chosen candidate models, with each model represented by a single parameter set $p$. A candidate model was generated by randomly assigning each parameter (${\displaystyle p_1,p_2,\stackrel{}{\dots }}$) a value comprised between 0% and 300% of its reference value ${\text{G}}_{\text{ref}}$. We then utilized a multivariate minimization algorithm (see Materials and methods for details) that evolves ${\displaystyle \mathbf{p}}$ until the GES optimization constraint defined by Equation 1 is satisfied. This method typically yields a large number of GES (7263 of the $\sim 10,000$ trials yield a GES with six parameters and three sensors described above, with 2737 either not converging or not producing a physiological output) and is computationally more efficient than a random search without optimization that yields very few GES.

Figure 3 with 3 supplements see all

Download asset Open asset

Results of the GES search are shown in Figure 3. Figure 3A shows the parameters of six representative GES and their corresponding AP waveforms (Figure 3B) and CaT profiles (Figure 3C). The CaT profiles are all very close to each other, which holds for all GES, while the AP waveforms exhibit larger variations owing to the fact that the GES search does not involve any voltage sensing. Figure 3D shows histograms of parameters for all GES. Conductances of sarcolemmal currents tend to be highly variable except for ${\displaystyle {\mathrm{I}}_{\mathrm{N}\mathrm{a}\mathrm{C}\mathrm{a}}}$, which turns out to be constrained by the intracellular sodium concentration sensor ($S_3={\left[\text{Na}\right]}_{\text{i}}$). This is revealed by a GES search with only Ca²⁺ sensing ($S_1$ and $S_2$) that yields a broader histogram for the ${\displaystyle {\mathrm{I}}_{\text{NaCa}}}$ conductance (Figure 3—figure supplement 1). The histogram of RyR expression level is very broad. This is consistent with the fact that this parameter was found to have a very weak effect on the CaT (Figure 2A,B) due to the compensatory adjustment of SR Ca²⁺ load (Figure 2E). In contrast, the histogram of SERCA is very narrow. This feature, which persists even if Na⁺ sensing is removed (Figure 3—figure supplement 1), is predominantly due to Ca²⁺ sensing. It stems from the fact that changing SERCA expression level has opposite effects on the CaT amplitude (Figure 2A) and average ${\left[\text{Ca}\right]}_{\text{i}}$ (Figure 2B), increasing one while decreasing the other or vice-versa. Therefore, those opposite effects cannot be compensated by changes of conductance of sarcolemmal currents that simultaneously increase or decrease both Ca²⁺ sensors, or by changes of RyR expression level that has a negligible effect on the CaT due to SR load adjustment. However, conductances of inward and outward currents that change both Ca²⁺ sensors in opposite directions can in principle compensate each other. This compensation is revealed by representing each GES as a point in a 3D plot (Figure 3E) whose axes are the conductances of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$ and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$. This plot shows that all GES lie close to a 2D surface in this 3D conductance space due to a three-way compensation between the effects of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$, and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ on the CaT. GES lie inside a smeared 2D surface (i.e. a 2D surface of finite thickness) in the 3D conductance space of Figure 3E. This feature stems from the fact that the GES parameter space considered here is in principle six-dimensional (four sarcolemmal current conductances and 2 Ca²⁺ protein expression levels). However, SERCA expression and ${\displaystyle {\mathrm{I}}_{\text{NaCa}}}$ conductance are constrained by Ca²⁺ and Na⁺ sensing, respectively, and RyR expression has a negligible effect on both Ca²⁺ and Na⁺ sensors, thereby reducing the relevant parameter space to the three conductance axes of Figure 3C. The subspace of GES that minimizes the cost function $E$ must therefore lie on the 2D surface $E=0$. This surface is smeared because ${\displaystyle {\mathrm{I}}_{\text{NaCa}}}$ is only constrained by Na⁺ sensing within a finite range and the GES search only minimizes $E$ within a finite tolerance ($E\le ϵ$ instead of $E=0$).

To facilitate the comparison with experiments presented in the next subsection, it is useful to represent the three-way compensation between ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$, and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ conductances by plotting the sum of the peak currents of ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$ and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ versus the peak current of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ with all three currents measured under voltage-clamp with a step from −50 to 0 mV. Those peak currents are proportional to conductances up to proportionality factors fixed by intra- and extracellular ionic concentrations and voltage. In this peak-current representation, the smeared 2D surface of GES of Figure 3E takes on the simpler form of a thick nearly straight line (Figure 3F). We note that even though correlations between two or more parameters have been explored in population models (Sánchez et al., 2014; Britton et al., 2013; Muszkiewicz et al., 2018), their sum studied here has not been previously considered.

In order to test the computational modeling predictions, and at the same time differentiate intra-heart cell-to-cell from inter-subject variability, we performed electrophysiological and contractile measurements on ventricular myocytes obtained from mouse hearts of nine different strains from the HMDP listed in the Materials and methods, each strain assumed to represent a different good enough solution. Peak values of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$, and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ were measured under voltage-clamp with a step from −50 to 0 mV following established protocols (see Materials and methods). The K⁺ currents were measured in the same cell and the Ca²⁺ currents in different cells. Contraction was analyzed by measuring mycoyte shortening during several paced beats in separate cells for six strains that include five of the strains in which conductances were measured. In order to collect enough statistics to distinguish cell-to-cell from inter-strain variability, several hearts of each isogenic strain were used. The number of cells that could be obtained from one heart for L-type Ca²⁺ current, several K⁺ currents, or contraction analysis varied from 1 to 7 so that several hearts of each strain were needed to obtain enough independent measurements to statistically distinguish intra-heart cell-to-cell from inter-strain variability (see data in Materials and methods).

The results of current and contraction measurements are shown in Figure 4. In Figure 4A, we plot the sum of the peak currents of ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$ and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ versus the peak current of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ together with the standard errors of the mean (SEM) of those quantities. Bar plots showing mean current values together with both SEM and standard deviation (SD) characterizing cell-to-cell variability are given in the Materials and methods. We also superimpose on this plot the computationally predicted GES of Figure 3F. Different HMDP strains, each representing a GES, are seen to function with different combinations of Ca²⁺ and K⁺ currents that compensate each other in a non-trivial three-way fashion that closely follows the GES computationally determined with a three-sensor search in which both the CaT and Na⁺ concentration are constrained (faded red points in Figure 4A). The sum of ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$ and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ follows a linear correlation with ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ (p=0.0007) using eight out of nine strains and the correlation remains statistically significant (p=0.0144) if the outlier strain (BXA12/PgnJ) is included. Interestingly, this outlier strain still falls within the larger range of computationally predicted GES using a two-sensor search without Na⁺ sensing (faded blue points in Figure 4A). To distinguish cell-to-cell from inter-strain variability, we performed a one-way ANOVA F-test on the ${\text{I}}_{\text{Ca},\text{L}}$ measurements. The result shows that ${\text{I}}_{\text{Ca},\text{L}}$ measurements for all strains do not originate from the same distribution (p-value p=0.000024). Furthermore, we performed a student T-test using raw data of ${\text{I}}_{\text{Ca},\text{L}}$ measurements for all pairs of strains. The results yield very small statistically significant p-values for pairs of strains with sufficiently different average current values (e.g. BXA25/PgnJ, CXB1/ByJ, and C57BL/6J in Figure 4A). Those results are consistent with the fact that mean currents differ much more than their standard error for those strains, as can be seen by visual inspection of means and SEM values corresponding to thin bars on both axes of Figure 4A. We conclude that inter-strain variability of ion channel conductances can be distinguished from cell-to-cell variability of those same quantities for a significant number of the strains investigated. While the Ca²⁺ and K⁺ currents were measured for the nine strains reported in Figure 4A, the Ca²⁺ current was measured in seven additional strains (total of 16 strains). Those additional measurements reported in the Materials and methods confirm that some strains can have markedly different ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ conductances.

Figure 4 with 1 supplement see all

Download asset Open asset

Unlike ion channel conductances, CaT properties were assumed not to vary in the computationally-enabled GES search, which rests on the hypothesis that Ca²⁺ sensing provides a feedback mechanism that regulates ion channel gene and protein expression. The results in Figure 4B, which use contraction as a surrogate for CaT amplitude, support this hypothesis by showing that mean values of cell shortening do not vary substantially across strains. This is confirmed by performing a standard ANOVA statistical test, which shows that cell shortening measurements for the six strains reported in Figure 4B do not originate from different distributions within statistical uncertainty (p-value p=0.4136).

Figure 5 with 2 supplements see all

Download asset Open asset

Current measurements discussed in the previous section show that mean conductances of Ca²⁺ and K⁺ currents vary between strains in a compensatory way so as to produce a normal CaT. They also show that conductances vary significantly from cell to cell around their mean values. This is illustrated in Figure 5A for three mouse strains that have statistically distinguishable mean ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ conductances (low, medium, and high) as measured by standard errors (thick bars), but exhibit large cell-to-cell variability as measured by the standard deviations (thin bars) of the distributions of conductance measurements in individual cells. This raises the question of whether compensation remains operative at the organ scale in the presence of large cell-to-cell variability. There are two interlinked aspects to this question. The first relates to the cellular-level dynamical coupling between membrane voltage and intracellular Ca²⁺ dynamics that is inherently nonlinear (Krogh-Madsen and Christini, 2012; Karma, 2013; Qu et al., 2014). Even when cells are uncoupled, this nonlinearity could potentially cause the mean CaT amplitude in an ensemble of cells with highly variable conductances to differ from the CaT amplitude computed in a single cell with conductances set to the mean values of the ensemble, as traditionally done in cardiac modeling. The second aspect relates to the additional effect of gap-junctional coupling between cells. This effect is well-known to smooth out cell-to-cell variation of AP waveforms on a mm scale that is much larger than the individual mycoyte length. However, whether this smoothing translates into increased organ-scale uniformity of CaT amplitude and contractility is unclear.

To address those two aspects, we constructed tissue scale computational models for three mouse strains with statistically distinguishable average conductances (as illustrated for ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ in Figure 5A). Tissues of each strain consisted of $56\times 56$ electrically coupled cells (see Materials and methods for details). Simulations were carried out with and without electrical coupling to assess the effect of the latter. The conductances of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$, and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ were assumed to vary randomly from cell to cell, and no constraint was imposed on the ratio of ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}+{\mathrm{I}}_{\text{Kur}}}$ to ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$, representing the worst case scenario in which both AP and CaT would exhibit maximal variations at the single myoycte level. Their values were drawn randomly from Gaussian distributions with average values and standard deviations that match experimental current measurements in each strain. All other parameters were kept fixed to reference values. The resulting cell-to-cell variation of conductances for three different strains is shown in Figure 5B) using the same peak-current representation of Figures 3F and 4A. In this representation, each point represents a different cell, and clouds of points of the same color represent all cells in a tissue of the same strain. Furthermore, the center of each cloud falls on the thick line corresponding to the computationally determined GES surface where compensation is operative at the single-cell level.

The results of simulations with populations of uncoupled and coupled cells with randomly varying conductances are shown in Figure 5C–G. Figure 5C shows that AP waveforms are highly variable when cells are uncoupled, reflecting the variability in conductances with no constraint imposed on the ratio of ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}+{\mathrm{I}}_{\text{Kur}}}$ to ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$. Figure 5D shows that AP waveforms becomes uniform when cells are coupled, as expected, even though interstrain variability is still significant. Figure 5E compares histograms of CaT amplitude and AP duration (APD) when cells are uncoupled and coupled. Consistent with the AP waveforms of Figure 5C and D, APD histograms in Figure 5E show that junctional coupling strongly reduces APD variability, as expected. CaT amplitude and average ${\displaystyle {\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$ histograms in turn reveal that, in coupled cells, the more uniform APD translates into a more uniform CaT amplitude and average ${\displaystyle {\left[\mathrm{C}\mathrm{a}\right]}_{\mathrm{i}}}$ (i.e. narrower $\Delta$Ca and <Ca> histograms, respectively), reflecting the influence of the cell’s APD on its CaT. At the organ scale, this ensures that cells in tissue have uniform APs. They also benefit modestly from a more uniform CaT as a result of coupling, promoting more uniform force generation throughout the tissue. Thus, tissue coupling compensates significantly when AP and CaT variability between single myoycytes is high, for the case in which the ratio of ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}+{\mathrm{I}}_{\text{Kur}}}$ to ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ is not constrained at the single myocyte level.

Figure 5F and G show that compensation remain operative at the organ scale. Figure 5F shows that, even though the strains have different average conductances (Figure 5B), they produce CaT amplitude histograms with approximately the same mean and width. In contrast, if the same simulation is repeated by fixing the conductance of K⁺ currents to the value of the strain with the medium value of ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ conductance (CXB1/ByJ), the different ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ conductances are not compensated by different ${\displaystyle {\mathrm{I}}_{\text{to},\text{f}}}$ and ${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ conductances, yielding CaT amplitude distributions with shifted peaks and hence different aggregate contraction (Figure 5G).

In summary, our results show that compensation remains operative at the organ scale because CaT amplitude histograms have similar means with and without electrical coupling (Figure 5E). This implies that cell-to-cell variability of conductances and hence APD causes variability of CaT amplitude without significantly affecting its mean, so that ${\displaystyle {\mathrm{I}}_{\text{Ca},\text{L}}}$ and potassium currents can compensate each other even though conductances exhibit large cell-to-cell variations from their mean values. Gap junctional coupling has the additional important effect of reducing CaT amplitude variability, thereby promoting tight organ-level behavior despite high cell-to-cell variability.

From a functional standpoint, the most relevant implication of the present study is that different GES may exhibit markedly different responses to perturbations, as previously demonstrated in a neuroscience context (Grashow et al., 2009). To examine this possibility, we reviewed data from separate studies of isoproterenol (ISO)-induced cardiac hypertrophy and heart failure in approximately 100 HMDP strains that include most of the strains used in the present study. In those studies, heart mass was measured in those strains before (${\text{m}}_{\text{pre}}$) and 3 weeks after (${\text{m}}_{\text{post}}$) implantation of a pump continuously delivering isoproterenol (Table 3). Figure 6 reveals the existence of a statistically very significant correlation between baseline ${\text{I}}_{\text{Ca},\text{L}}$ conductance and the hypertrophic response (${\text{m}}_{\text{post}}$/${\text{m}}_{\text{pre}}$). Although many factors contribute to the hypertrophic response in the HMDP (Rau et al., 2017; Santolini et al., 2018), intracellular Ca²⁺ overload activating the Ca²⁺-calcineurin-NFAT signaling pathway has been shown to play a major role (Bers, 2008). Since ${\text{I}}_{\text{Ca},\text{L}}$ is the major pathway of Ca²⁺ entry into the cytoplasm, it is intriguing to speculate that strains with a larger baseline ${\text{I}}_{\text{Ca},\text{L}}$ conductance under baseline conditions have a more robust increase in ${\text{I}}_{\text{Ca},\text{L}}$ that is not adequately compensated by repolarizing ${\text{K}}^{+}$ currents, making those strains more susceptible to Ca²⁺ overload when ${\text{I}}_{\text{Ca},\text{L}}$ is enhanced during sustained $\beta$-adrenergic stimulation by isoproterenol. Hypothetically, this may result in a stronger cardiac hypertrophic response. To make this case convincingly, however, would require demonstrating that Ca²⁺ overload is chronically worsened in strains with a high baseline ${\text{I}}_{\text{Ca},\text{L}}$ conductance and ruling out other strain-dependent hypertrophy-promoting pathways that are not Ca²⁺-dependent, which is beyond the scope of the present work.

Figure 6

Download asset Open asset

In a neuroscience context, ionic conductances of neurons from the stomatogastric ganglion of different crabs were previously found by Schulz et al. (2006) to be correlated with gene expression, as shown by independent measurements in the same subjects of functional densities of different ion channels, used to determine conductances, and mRNA levels of genes encoding for pore-forming subunits of those channels. In a cardiac context, decrease of ${\text{I}}_{\text{Ca},\text{L}}$ current density has been shown to be correlated with a decrease of Cav1.2 mRNA expression in response to a sustained increase of pacing rate in cultured adult canine atrial cardiomyocytes mimicking atrial tachycardia remodeling (Qi et al., 2008). In the present study, we did not perform independent measurements of gene expression in the same ventricular myocytes used to measure ionic conductances. However, to examine the possible relationship between compensation of conductances and gene expression, we reviewed the gene expression data from the aforementioned studies of ISO-induced cardiac hypertrophy and heart failure in approximately 100 HMDP strains that include most of the strains used in the present study. Gene expression was measured both in control (pre-ISO) and after injection of ISO for 21 days in 8- to 10-week-old female mice (post-ISO). Details of heart biopsies conducted pre- and post-ISO and microarray data analysis are given in the Methods section of Santolini et al. (2018).

No statistically significant pairwise correlation (Pearson correlation coefficient $r<0.25$ and p-value $p>0.05$) were found between expression levels of the genes Cacna1c, Kcnd2, and Kcna5, encoding for the pore forming subunit of the Cav1.2, Kv4.2, and Kv1.5 channels associated with ${\text{I}}_{\text{Ca},\text{L}}$, ${\text{I}}_{\text{to},\text{f}}$, and ${\text{I}}_{\text{Kur}}$, respectively. However, a statistically very significant correlation was found between expression levels of Cacna1c and Kcnip2 that encodes the KChIP2 accessory $\beta$ subunits directly interacting with Kv4.2 (see Figure 7 and caption for $r$ and $p$ values). This correlation is present both in control (pre-ISO), which is the condition relevant to our conductance measurements in selected strains (Figure 4), and post-ISO. Since increased KChiP2 level is known to increase the functional current density of ${\text{I}}_{\text{to},\text{f}}$ (Kuo et al., 2001; Jin et al., 2010), the strong positive correlation between Cacna1c (Cav1.2) and Kcnip2 (KChIP2) expression levels may partially contribute to the positive correlation between ${\text{I}}_{\text{Ca},\text{L}}$ and ${\text{I}}_{\text{to},\text{f}}$+${\displaystyle {\mathrm{I}}_{\text{Kur}}}$ functional current densities (Figure 4).

Figure 7

Download asset Open asset

The hybrid mouse diversity panel (HMDP) consists of a population of over 100 inbred mouse strains selected for usage in systematic genetic analyses of complex traits (Ghazalpour et al., 2012). The main goals in selecting the strains were to (i) increaseresolution of genetic mapping, (ii) have a renewable resource that is available to all investigators world-wide, and (iii) provide a shared data repository (https://systems.genetics.ucla.edu/about/hmdp2) that would allow the integration of data across multiple scales, including genomic, transcriptomic, metabolomic, proteomic, and clinical phenotypes.

At sacrifice, hearts were excised, drained of excess blood and weighed. Each chamber of the heart (LV with inter-ventricular septum, RV-free wall, RA and LA) was isolated and subsequently weighed. Cardiac hypertrophy was calculated as the increase in total heart weight after isoproterenol (ISO) treatment compared to control animals (see Table 3). As described prevoiously (Wang et al., 2016), the ISO treatment consisted of 30 mg per kg body weight per day of Isoproterenol (ISO) administered for 21 days in 8- to 10-week-old female mice using ALZET osmotic mini-pumps, which were surgically implanted intraperitoneally. The average number of control hearts per strain was 2.75. The average number of treated hearts per strain was 3.5. The exact number of hearts per strain can be found in Table S1 of Santolini et al. (2018).

We have developed a novel mathematical model of mouse ventricular myocytes that combines elements of previously published ventricular mycoyte models (Shiferaw et al., 2003; Shannon et al., 2004; Bondarenko et al., 2004; Mahajan et al., 2008). For this purpose, we kept the mathematical formulation of intracellular calcium cycling of the Mahajan model developed by Shiferaw et al. (2003), which physiologically incorporates graded release by linking the Ca²⁺ spark recruitment rate to ${\displaystyle {\text{I}}_{\text{Ca},\text{L}}}$ current magnitude, and replaced several sarcolemmal currents by those formulated by Bondarenko et al., 2004 for mouse ventricular mycoytes. This model allowed us to explore efficiently the space of good enough solutions that we can compare to experimentally measured variability found in the hybrid mouse diversity panel (HMDP). The Bondarenko et al., 2004 mouse model has sarcolemmal currents fitted to detailed experimental measurements of sarcolemmal currents in mouse ventricular myocytes. The Mahajan et al. (2008) model is a rabbit model that integrates a Markov model of ${\displaystyle {\text{I}}_{\text{Ca},\text{L}}}$ together with the Shannon et al. (2004) formulation of other sarcolemmal currents and the Shiferaw et al. (2003) model of calcium cycling and SR calcium release. The Shiferaw et al. (2003) model represents the release of calcium from the SR as a sum of individual spark events, which reproduces important observed instabilities such as Ca²⁺ transient alternans. Even though the Shiferaw et al. (2003) model of calcium cycling model was developed for rabbit mycoytes, it can in principle describe Ca²⁺ cycling in other species including mouse. Therefore, we constructed a mouse model starting with the Mahajan et al. (2008) model, which incorporates the Shiferaw et al. (2003) model of calcium cycling, and using the formulation of Shannon et al. (2004) for ${\displaystyle {\text{I}}_{\text{Ca},\text{L}}}$ and Bondarenko et al., 2004 for other sarcolemmal currents fitted to mouse data. We did not use the Bondarenko et al., 2004 detailed Markov formulations of the gating of ${\displaystyle {\text{I}}_{\text{Na}}}$ and ${\displaystyle {\text{I}}_{\text{Ca},\text{L}}}$ channels with fast transition rates that were computationally prohibitive for the number of simulations we performed in this study. However, we checked that our combined model reproduces well the mouse electrophysiological phenotype of the Bondarenko et al., 2004 model while being computationally efficient and incorporating a realistic description of Ca²⁺ cycling. We verified that the Shiferaw et al. (2003) model of Ca²⁺ cycling integrated with the Shannon et al. (2004) model of ${\text{I}}_{\text{Ca},\text{L}}$ produced a normal bell-shaped SR release as a function of step-voltage. We also verified that the force-frequency relationship produced by the model has the correct negative staircase observed experimentally. The equations for the model are described below. The reference values of the six parameters varied in this study are given in Table 4 and produce baseline AP and CaT morphologies consistent with the experimental measurements in Bondarenko et al., 2004. Table 5 lists all other parameters used that were kept fixed.