The Ca2+ transient as a feedback sensor controlling cardiomyocyte ionic conductances in mouse populations

Conductances of ion channels and transporters controlling cardiac excitation may vary in a population of subjects with different cardiac gene expression patterns. However, the amount of variability and its origin are not quantitatively known. We propose a new conceptual approach to predict this variability that consists of finding combinations of conductances generating a normal intracellular Ca2+ transient without any constraint on the action potential. Furthermore, we validate experimentally its predictions using the Hybrid Mouse Diversity Panel, a model system of genetically diverse mouse strains that allows us to quantify inter-subject versus intra-subject variability. The method predicts that conductances of inward Ca2+ and outward K+ currents compensate each other to generate a normal Ca2+ transient in good quantitative agreement with current measurements in ventricular myocytes from hearts of different isogenic strains. Our results suggest that a feedback mechanism sensing the aggregate Ca2+ transient of the heart suffices to regulate ionic conductances.


Introduction
Following the landmark publication of the Hodgkin-Huxley model of nerve-cell action potential over six decades ago (Hodgkin and Huxley, 1952), electrophysiological models of increasing complexity have been developed to describe the cardiac action potential (AP) and its interaction with the intracellular calcium (Ca 2+ ) signal (Noble, 2011;Silva and Rudy, 2010), which links electrical signaling to mechanical contraction in cardiomycoytes (Bers, 2001). As illustrated in Figure 1 for a mouse ventricular mycoyte (Bondarenko et al., 2004), those models typically involve a large set of interacting cellular components that includes various voltage-gated membrane ion channels and transporters, the Na + /Ca 2+ exchanger, and Ca 2+ handling proteins such as the ryanodine receptor (RyR) Ca 2+ release channels, which open in response to Ca 2+ entry into the cell via L-type Ca 2+ channels, and the sarcoplasmic reticulum (SR) Ca 2+ ATPase (SERCA), which uptakes Ca 2+ back into the SR. Ca 2+ release and uptake from the SR causes a transient rise in cytosolic Ca 2+ concentration, the calcium transient (CaT), which activates myocyte contraction. Those cellular-scale models have been traditionally constructed by piecing together separate mathematical models describing molecular-scale components in the same species (guinea pig (Luo and Rudy, 1991), mouse (Bondarenko et al., 2004), rabbit (Shannon et al., 2004;Mahajan et al., 2008), dog (Fox et al., 2002), etc.), and by sometimes mixing models in different species. Experimental measurements of voltage-current relationships and other properties used to develop those models are typically averaged over several cells from one or a few hearts. While such prototypical 'population-averaged' (and even 'speciesaveraged') models have proven useful to investigate basic mechanisms of cardiac arrhythmias on cellular to tissue scales (Krogh-Madsen and Christini, 2012;Karma, 2013;Qu et al., 2014), they fall short of predicting how different individuals in a genetically diverse population respond to perturbations (such as physiological stressors, ion channel mutations, drug or gene therapies, etc.) affecting one or several cellular components.
In a neuroscience context, the limitation of population-averaged models has been highlighted by pioneering theoretical and experimental studies by Abbot, Marder, and co-workers demonstrating that ion channel conductances can exhibit a high degree of activity-dependent plasticity as well as variability between individuals (LeMasson et al., 1993;Siegel et al., 1994;Liu et al., 1998;Golowasch et al., 1999;Prinz et al., 2004;Schulz et al., 2006;Marder and Goaillard, 2006;Grashow et al., 2009;Marder, 2011;O'Leary et al., 2013). Theoretical studies along this line first originated from an attempt to explain why neurons can maintain fixed electrical activity patterns despite a high rate of ion channel turnover. This question was addressed by treating ion channel conductances as dynamical variables in models of neuronal activity and by using the intracellular Ca 2+ concentration as an activity-dependent feedback mechanism regulating their values (LeMasson et al., 1993;Siegel et al., 1994;Liu et al., 1998), a mechanism supported by experiments (Golowasch et al., 1999). Those model neurons displayed remarkable properties such as the ability to modify their conductances to maintain a given behavior when perturbed or to develop different properties in response to different patterns of presynaptic activity. Subsequently, a different type of computational study in which model parameters (conductances and synaptic strengths of a circuit model of the crustacean somato-gastric ganglion) were varied randomly, demonstrated that a similar bursting activity could be obtained with multiple parameter sets, that is multiple 'good enough solutions' (GES) (Prinz et al., 2004). This prediction agreed qualitatively with an experimental study demonstrating that neurons of the same circuits obtained from different crabs could have markedly different ion channel densities, corresponding to different gene expression, but yet different circuits could generate similar bursting activities (Schulz et al., 2006).
A later experimental study further showed that different circuits, while operating similarly under controlled conditions, could respond differently to perturbations such as serotonin addition, which increases the bursting frequency at a population level, but lowers it in some individuals (Grashow et al., 2009). Those findings shed light on why pharmacological treatments may work in some individuals but not others. In addition, they suggest that the existence of different good enough solutions provides an evolutionary advantage for the survival of a genetically diverse population by allowing different individuals to better adapt to different environmental challenges so as to survive and restore the population.
Those studies and related ones in the broader context of systems biology (Daniels et al., 2008;Gutenkunst et al., 2007;Transtrum et al., 2015) and cardiac electrophysiology (Sarkar and Sobie, 2009;Sarkar and Sobie, 2010;Weiss et al., 2012;Sarkar et al., 2012;Britton et al., 2013;Groenendaal et al., 2015;Muszkiewicz et al., 2016;Krogh-Madsen et al., 2016) have produced a paradigm shift from population-averaged models, with unique fine-tuned parameter sets, to populations of models characterized by multiple parameter sets. In this new paradigm, each set representing a different individual in a population can produce a similar behavior under controlled conditions, but a starkly different response to perturbations for some individuals. This paradigm shift, however, creates new theoretical and experimental challenges.
On the theoretical side, an open question is how to search for GES representing different individuals in a population. The results of this search, which has been carried out using various methods (e.g. random search (Prinz et al., 2004), multivariate regression analysis (Sarkar and Sobie, 2009;Sarkar and Sobie, 2010), or genetic algorithms (Groenendaal et al., 2015), depend critically on what outputs are selected to constrain model parameters. To date, parameter searches in a cardiac context have been 'AP centric', focusing primarily on features of the membrane voltage (V m ) signal. Sarkar and Sobie showed that very different combinations of ion conductances can produce almost identical cardiac AP waveforms, albeit different CaT amplitudes, and that adding additional constraints on the V m and Ca 2+ signals can further constrain model parameters (Sarkar and Sobie, 2009;Sarkar and Sobie, 2010). Groenendaal et al., 2015 constrained model parameters using V m traces with variable AP waveforms recorded from guinea pig cardiomyocytes under randomly timed electrical stimuli, as opposed to a unique AP waveform recorded during periodic pacing. This search yielded parameter sets that are potentially better suited to describe more complex aperiodic forms of V m dynamics relevant for arrhythmias. Britton et al., 2013 observed experimentally a significant variability in AP waveform in rabbit Purkinje fibers and searched for model parameter combinations that reproduce this waveform variability. They then used those parameter sets to predict different effects of pharmacological blockade of cardiac HERG (I Kr current) potassium channel in different subjects (Britton et al., 2013).
All those GES searches have relied for the most part on using measured AP features (Sarkar and Sobie, 2009;Sarkar and Sobie, 2010;Sarkar et al., 2012;Britton et al., 2013;Groenendaal et al., 2015;Muszkiewicz et al., 2016;Krogh-Madsen et al., 2016) to constrain model parameters, even though some recent studies have also considered the additional effect of constraining the CaT Mayourian et al., 2017). However, given that there is no known voltage-sensing mechanisms regulating ion channel expression, it remains unclear if natural biological variability can be predicted based on AP features. Here, we adopt a different 'Ca 2+ centric' view (Weiss et al., 2012), which postulates as in a neuroscience context (Golowasch et al., 1999;LeMasson et al., 1993;Siegel et al., 1994;Liu et al., 1998;O'Leary et al., 2013) that model parameters are predominantly constrained by feedback sensing of Ca 2+ , and potentially other ions (e.g. Na + ) affecting ion channel regulation. Our hypothesis is that the CaT is critical for generating blood pressure, which is sensed by the carotid baroreceptors and feeds back through the autonomic nervous system to regulate the CaT via controlling levels of Ca-cyling proteins and the AP in a way that preserves blood pressure. This provides a very straightforward physiological mechanism that we show not only constrains the CaT to a physiological waveform, but, as an added and novel bonus, also constrains AP features through the ratio of inward Ca currents and outward K currents. Under this hypothesis, multiple parameter combinations producing a normal CaT could potentially represent different GES in a genetically diverse population. In addition, unlike voltage, intracellular concentrations of Ca 2+ and Na + ions ( Ca ½ i and Na ½ i , respectively) have a known interactive role in transcriptional regulation of cardiac ion channel proteins and their function (Rosati and McKinnon, 2004). For example, the Ca 2+ /calcineurin/NFAT pathway regulates L-type Ca 2+ channel (LCC) expression (Qi et al., 2008) and Na + modulates cAMP-dependent regulation of ion channels in the heart (Harvey et al., 1991) including phosphorylation of LCCs via cAMP-dependent protein kinase (Balke and Wier, 1992). To test this hypothesis, we perform a GES search in which parameters of a mouse ventricular myocyte model are only constrained by CaT features and Na ½ i . This search yields GES with different conductances of the L-type Ca 2+ current (I Ca;L ) and K + currents (I to;f and I Kur ) and reveals that conductances are strongly correlated due to compensatory effects of those currents on the CaT.
On the experimental side, a major challenge is to test whether different GES produced by any given search method are representative of different individuals in a genetically diverse population. Performing this test generally requires distinguishing quantitatively the variability of conductances and electrophysiological phenotype observed in cells extracted from the same heart (intra-heart cellto-cell variability) from the variability of the same quantities between different subjects (inter-subject variability). Making this distinction is made extremely difficult by the fact that AP features and conductances vary significantly between cells extracted from same region of the heart (Banyasz et al., 2011;Groenendaal et al., 2015) and that regional (e.g. ventricular base-to-apex and epicardium to endocardium) variations of ion channel expression are also present. The existence of large intra-heart cell-to-cell variability, and the practical limitation that only a finite number of cells can typically be extracted from a single heart for current measurements, raises the question of whether it is feasible to distinguish electrophysiological parameters between genetically distinct individuals.
To cope with this challenge, we use here the Hybrid Mouse Diversity Panel (HMDP) that is a collection of approximately 100 well-characterized inbred strains of mice that can be used to analyze the genetic and environmental factors underlying complex traits. Because inbred strains are isogenic and renewable, we are able to use multiple hearts from the same strain to obtain enough statistics to differentiate quantitatively between intra-heart and inter-subject variability in conductances of key currents (I Ca;L , I to;f and I Kur ) affecting the AP and CaT of mouse ventricular myocytes from different strains. The results show that, despite large cell-to-cell variability, some strains have clearly distinguishable mean conductances (i.e. conductances averaged over all cells for the same strain). Mean conductances can vary by up to two-and-a-half fold between strains. The results further show that, remarkably, variations of mean conductances for individual strains follow the same correlation (I Ca;L current is large or small when the sum of I to;f and I Kur currents is large or small, respectively) predicted by our computational Ca 2+ centric GES search. The central hypothesis that parameters are constrained predominantly by features of the CaT (as a surrogate for arterial blood pressure) is further validated experimentally by showing that strains with very different conductances have similar contractile activity. It is worth emphasizing that the main novelty of the present study is the use of the HMDP to validate this hypothesis. The computational identification of GES itself uses a standard search algorithm, which consists of minimizing a cost function constructed from features of the CaT and the intracellular sodium concentration. In addition, we use tissue scale simulations to show that compensation remains effective at an organ scale despite large cell-to-cell variability within an individual heart. Finally, we use our findings to interpret the results of recent studies of cardiac hypertrophy and heart failure induced by a stressor in the HMDP (Ghazalpour et al., 2012;Rau et al., 2015;Rau et al., 2017;Santolini et al., 2018).

Effects of individual conductances on the calcium transient
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 D Ca ½ i between the diastolic and peak value of the cytosolic Ca 2+ concentration Ca ½ i , and the time-averaged value of Ca ½ i over one pacing period, denoted by h Ca ½ i i. The CaT amplitude D Ca ½ i is a major determinant of the contractile force while h Ca ½ i i provides an average measure of the cytosolic Ca 2+ concentration in the cell. Both quantities will be used as Ca 2+ 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 2+ handling proteins that include the ryanodine receptor (RyR) Ca 2+ release channels and the sarcoplasmic reticulum (SR) Ca 2+ ATPase SERCA, which pumps Ca 2+ 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 D Ca ½ i and h Ca ½ i i and the intracellular sodium concentration Na ½ i no longer vary from beat to beat. Figure 2 shows the effects of individual parameter changes on the steady-state CaT amplitude ( Figure 2A) and average Ca ½ 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 2+ current I Ca;L is seen to strongly increase both D Ca ½ i and h Ca ½ i i but has a  weak effect on the AP waveform ( Figure 2C). The effect on the CaT stem from the fact that I Ca;L is the main trigger of Ca 2+ -induced Ca 2+ release(CICR), which transfers a large amount of Ca 2+ from the SR to the cytosol. The weak effect on APD is due to the fact that the increase of CaT amplitude causes I Ca;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 I to;f (fast inactivating component of the transient outward current) and I Kur causes both D Ca ½ i and h Ca ½ i i to decrease. Increasing either of those K + currents speeds up repolarization (as illustrated for I to;f in Figure 2D) and hence inactivation of I Ca;L , thereby reducing the magnitude of CICR. Increasing the conductance of the sodium-calcium exchanger current I NaCa also causes both D½Ca i and h Ca ½ i i to decrease by enhancing the forward mode of this current that extrudes Ca 2+ 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 I Ca;L conductance on CaT amplitude in the present mouse model. Changing RyR expression from its reference value is seen to leave D Ca ½ i and h Ca ½ i i almost unchanged, even though it strongly affects the SR Ca 2+ concentration Ca ½ 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 2+ release flux through RyRs) yields a transient increase in CaT amplitude, but no change in the steady-state CaT amplitude after Ca ½ SR adjusts to a lower steady-state level (Bers, 2001). This effect is illustrated by time traces of Ca ½ SR and Ca ½ i in steady-state for different RyR expression levels in Figure 2E. Finally, changing the expression level of SERCA has opposite effects on D Ca ½ i and h Ca ½ i i. Increasing SERCA magnitude increases SR Ca 2+ load, thereby increasing the amount of SR Ca 2+ release and CaT amplitude, but at the same time depletes Ca 2+ from the cytosol.

Computationally determined good enough solutions
Next, we performed a computational search for combinations of parameters that yield a normal electrophysiological output as defined by the steady-state CaT amplitude D Ca ½ i , time averaged cytosolic Ca 2+ concentration h Ca ½ i i, and intracellular sodium concentration Na ½ i at a 4 Hz pacing frequency. A GES search that uses the diastolic and peak Ca ½ i values as Ca 2+ sensors, instead of D Ca ½ i (the difference between the peak and diastolic Ca ½ i values) and time averaged Ca ½ i , gives nearly identical results. So our Ca 2+ 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 D Ca ½ * i , h Ca ½ i i Ã , and Na ½ * i 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 p ð Þ from their desired target values S * n . Here, and is a small tolerance that we choose to be 5%. E is a function of model parameters p ¼ p 1 ; p 2 ; . . . ð Þ chosen to consist of the conductances of I Ca;L , I to;f , I Kur , and I 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 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 (p 1 ; p 2 ; . . .) a value comprised between 0% and 300% of its reference value G ref . We then utilized a multivariate minimization algorithm (see Materials and methods for details) that evolves p until the GES optimization constraint defined by Equation 1 is satisfied. This method typically yields a large number of GES (7263 of the~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.
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 I NaCa , which turns out to be constrained by the intracellular sodium concentration sensor (S 3 ¼ Na ½ i ). This is revealed by a GES search with only Ca 2+ sensing (S 1 and S 2 ) that yields a broader histogram for the I 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 2+ 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 2+ sensing. It stems from the fact that changing SERCA expression level has opposite effects on the CaT amplitude ( Figure 2A) and average Ca ½ 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 2+ 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 2+ 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 I Ca;L , I to;f and I Kur . This plot shows that all GES lie close to a 2D surface in this 3D conductance space due to a threeway compensation between the effects of I Ca;L , I to;f , and I 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 2+ protein expression levels). However, SERCA expression and I NaCa conductance are constrained by Ca 2+ and Na + sensing, respectively, and RyR expression has a negligible effect on both Ca 2+ 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 I NaCa is only constrained by Na + sensing within a finite range and the GES search only minimizes E within a finite tolerance (E 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 I Ca;L , I to;f , and I Kur conductances by plotting the sum of the peak currents of I to;f and I Kur versus the peak current of I Ca;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.

Good enough solutions in the HMDP
In order to test the computational modeling predictions, and at the same time differentiate intraheart 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 I Ca;L , I to;f , and I 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 2+ 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 2+ 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 I to;f and I Kur versus the peak current of I Ca;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 2+ 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 I to;f and I Kur follows a linear correlation with I Ca;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 I Ca;L measurements. The result shows that I Ca;L measurements for all strains do not originate from the same distribution (pvalue p=0.000024). Furthermore, we performed a student T-test using raw data of I Ca;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 2+ and K + currents were measured for the nine strains reported in Figure 4A, the Ca 2+ 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 I Ca;L conductances.
Unlike ion channel conductances, CaT properties were assumed not to vary in the computationally-enabled GES search, which rests on the hypothesis that Ca 2+ 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). between conductances of I Ca;L , I to;f , and I Kur . Each GES is represented by a red dot. All GES lie close to a 2D surface in this 3D plot. Pairwise projections (grey shadows) do not show evidence of two-way compensation between pairs of conductances. (F) Alternate representation of three-way compensation obtained by plotting I Ca;L versus the sum of I to;f and I Kur . Peak values of those currents after a voltage step from À50 to 0 mV are used to make this plot that can be readily compared to experiment. Different time windows are plotted for the AP waveforms and CaT in B and C, respectively. DOI: https://doi.org/10.7554/eLife.36717.004 The following figure supplements are available for figure 3: showing quantitative agreement between theoretically predicted and experimentally measured compensation of inward Ca 2+ and outward K + currents. Equivalent plot of Figure 3F showing the sum of I to;f and I Kur versus I Ca;L for nine different mouse strains using peak values of those currents (proportional to conductances) after a voltage step from À50 to 0 mV. Mean current values (green filled squares) are shown together with standard errors of the mean (thin bars) for each strain. The number of cells used for each strain is given in Table 1 of the Materials and methods section. Computationally determined GES are superimposed and shown as faded red points using all three sensors (CaT amplitude, average Ca ½ i , and diastolic Na ½ i ) and faded blue points for two sensors (CaT amplitude and average Ca ½ i ). Lines represent linear regression fits using the method of Chi-squared minimization with errors in both coordinates including (solid line, p=0.0144) and excluding (dashed line, p=0.0007) the outlier strain BXA12/PgnJ marked by a red box. The small p values of those fit validate the computationally predicted three-way compensation of Ca 2+ and K + currents. The three strains selected for the organ scale study (C57BL/6J, CXB1/ByJ, and BXA25/PgnJ) with low, medium, and high I Ca;L conductance, respectively, are highlighted by blue circles. (B) Cell shortening, measured as the fraction of resting cell length at 4 Hz pacing frequency in different HMPD strains where thick and thin bars correspond to standard error of the mean and standard deviation, respectively. A standard ANOVA test shows no significant differences in cell shortening between strains (p=0.4136) supporting the hypothesis that different combinations of conductances produce a similar CaT and contractile activity. Figure 4 continued on next page Compensation at the organ scale Current measurements discussed in the previous section show that mean conductances of Ca 2+ 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 I Ca;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 2+ 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 I Ca;L in Figure 5A). Tissues of each strain consisted of 56 Â 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 I Ca;L , I to;f , and I Kur were assumed to vary randomly from cell to cell, and no constraint was imposed on the ratio of I to;f þ I Kur to I Ca;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 I to;f þ I Kur to I Ca;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 Ca ½ i histograms in turn reveal that, in coupled cells, the more uniform APD translates into a more uniform CaT amplitude and average Ca ½ i (i.e. narrower DCa 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 I to;f þ I Kur to I Ca;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 I Ca;L conductance (CXB1/ByJ), the different I Ca;L conductances are not compensated by different I to;f and I 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 I Ca;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.
Cardiac hypertrophic response to a stressor 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 (m pre ) and 3 weeks after (m post ) implantation of a pump continuously delivering isoproterenol ( Table 3). Figure 6 reveals the existence of a statistically very significant correlation between baseline I Ca;L conductance and the hypertrophic response (m post /m pre ). Although many factors contribute to the hypertrophic response in the HMDP (Rau et al., 2017;Santolini et al., 2018), intracellular Ca 2+ overload activating the Ca 2+ -calcineurin-NFAT signaling pathway has been shown to play a major role (Bers, 2008). Since I Ca;L is the major pathway of Ca 2+ entry into the cytoplasm, it is intriguing to speculate that strains with a larger baseline I Ca;L conductance under baseline conditions have a more robust increase in I Ca;L that is not adequately compensated by repolarizing K þ currents, making those strains more susceptible to Ca 2+ overload when I Ca;L is enhanced during sustained b-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 2+ overload is chronically worsened in strains Figure 5 continued I Ca;L conductance, respectively, and the grey points are the results of the three-sensor GES search (same as Figure 3F). (C) Variable AP waveforms in uncoupled myocytes with conductances randomly chosen from the distribution shown in B for C57BL/6J and D) AP waveforms for coupled myocytes in tissue for C57BL/6J and the two other strains. AP waveforms of uncoupled cells vary significantly from cell to cell as observed experimentally (Fig.  Figure 5-figure supplement 1) but are uniform in electrotonically coupled cells, as expected. (E) Histograms of Ca 2+ transient (CaT) amplitude (DCa) and action potential duration (APD) for C57BL/6J in electrotonically uncoupled and coupled cells. Importantly, in coupled cells, the more uniform APD translates into a much more uniform CaT amplitude, reflecting the strong effect of the cell's APD on its CaT amplitude. (F) Distribution of CaT amplitudes within electrotonically coupled cells in tissue scale simulations using the parameter distributions from B. The three strains have the same mean CaT amplitude averaged over all cells marked by a thick vertical gray line, thereby demonstrating that compensation of Ca 2+ and K + currents remains operative at a tissue scale. (G) Distribution of CaT amplitudes obtained by varying only I Ca;L conductance and with I to;f and I Kur conductances fixed to their reference values. Lack of compensation between Ca 2+ and K + currents in this case yields different mean CaT amplitude. DOI: https://doi.org/10.7554/eLife.36717.011 The following figure supplements are available for figure 5: with a high baseline I Ca;L conductance and ruling out other strain-dependent hypertrophy-promoting pathways that are not Ca 2+ -dependent, which is beyond the scope of the present work.

Compensation and gene expression
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 I Ca;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 I Ca;L , I to;f , and I Kur , respectively. However, a statistically very significant correlation was found between expression levels of Cacna1c and Kcnip2 that encodes the KChIP2 accessory b 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 I to;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 I Ca;L and I to;f +I Kur functional current densities (Figure 4).

Discussion
In the present study, we have proposed a new methodology for searching for combinations of electrophysiological parameters representing different individuals in a genetically diverse population. While previous studies have used primarily AP features to constrain parameters (Sarkar and Sobie, 2009;Sarkar and Sobie, 2010;Sarkar et al., 2012;Britton et al., 2013;Groenendaal et al., 2015;Muszkiewicz et al., 2016;Krogh-Madsen et al., 2016), we have chosen to constrain parameters using the Ca 2+ transient that plays a key role to regulate ion channel expression and activity. This choice is based on a straightforward physiological hypothesis, namely that the CaT is critical for generating blood pressure, which is sensed by the carotid baroreceptors and feeds back through the autonomic nervous system to regulate the CaT in a way that preserves blood pressure. In contrast, a physiological basis for sensing cardiac voltage to regulate the AP and CaT is unclear. We have also examined the effect of additionally constraining the intracellular Na + concentration that is also known to modulate ion channel activity. Regulatory mechanisms traverse different levels of biological organization from transcriptional regulation to post-transcriptional and posttranslational modification to ion channel trafficking and phosphorylation. Those mechanisms are presently not known in sufficient detail to be modeled quantitatively. However, there is sufficient experimental evidence of feedback sensing of cellular activity via Ca 2+ (Qi et al., 2008) and Na + (Harvey et al., 1991;Balke and Wier, 1992) concentrations to make a search that constrains model parameters based on those signals plausible. The Ca 2+ transient determines the contractile force underlying arterial blood pressure generation regulated by baroreceptor feedback via the autonomic nervous system. Hence, fixing the diastolic and peak Ca ½ i values is a physiologically meaningful choice to search for parameter combinations that produce a normal diastolic and systolic function, which we have adopted here. Previous work  has provided evidence of a compensatory increase of I Ks following exposure of canine cardiomyocytes to a pharmacological I Kr blocker. Even though the mechanisms are not clear, it has been hypothesized that feedback sensing of Ca ½ i may potentially underlie the compensatory upregulation of I Ks through post-transcriptional upregulation of underlying channel subunits mediated by microRNA changes. Together with other studies (Qi et al., 2008), those previous findings may provide supportive evidence for the present Ca 2+ sensing hypothesis and suggests its generality beyond mice.
A remarkable and nontrivial finding of the present computational study is that Ca 2+ sensing suffices to produce a physiological AP waveform whose duration spans comparable range (see Figure 5C) to that recorded experimentally in isolated myocytes ( Figure 5-figure supplement 1), even though the voltage signal is not used to constrain model parameters. For comparison, we show in Figure 3-figure supplement 3 the results of a GES search that uses voltage instead of Ca 2+ sensing. With voltage sensing alone (both without and with Na ½ i sensing), the AP waveform was readily constrained as expected, but the CaT was highly variable and often not physiological. Moreover, the correlation between inward Ca 2+ and outward K + currents observed in the HMDP ( Figure 4A) was no longer preserved, since other inward and outward currents including I NaCa could regulate AP duration when the CaT was not constrained. Those results support our hypothesis that ionic conductances are primarily regulated by feedback mechanisms sensing ionic concentrations. Since the CaT is critical for generating blood pressure (which is sensed by the carotid baroreceptors and feeds back through the autonomic nervous system to regulate the CaT by controlling levels of Ca-cyling proteins and the AP in a way that preserves blood pressure), Ca 2+ sensing provides a straightforward physiological mechanism that not only constrains the CaT to a physiological waveform, but also, as an added and novel bonus, constrains AP features through the ratio of inward Ca 2+ current and outward K + currents. It is much less clear, on the other hand, how voltage would be sensed by the heart to provide a feedback mechanism to control both the AP waveform and the CaT. Voltage-sensing alone provided a reliable AP waveform, but a highly unreliable CaT as shown in Figure 3-figure supplement 3.
Our results show that CaT calibration results in considerable AP waveform variability ( Figure 5C) in isolated myocytes as expected if the AP waveform is not constrained. This finding is consistent with the observations that AP variability is considerable when measured experimentally in patch clamp studies ( Figure 5-figure supplement 1), but greatly reduced in tissue because less frequent atypical AP waveforms are voltage-clamped by the more typical AP waveforms of their neighbors.
Our finding that CaT calibration results in considerable AP waveform variability is also consistent with the converse finding in a previous study (Muszkiewicz et al., 2018) and here (Figure 3-figure  supplement 3) that the CaT is highly variable when the AP waveform alone is constrained without Figure 7. Compensation and gene expression. Plot showing the existence of a statistically very significant correlation (Pearson correlation coefficient r ¼ 0:47 and p-value, p ¼ 8:1 10 À13 ) between the expression level of Kcnip2, encoding the KChIP2 accessory b subunits that interact with Kv4.2 channels (I to;f ) and of Cacna1c, a gene encoding the a1C subunit of the Cav1.2 L-type calcium channels (I Ca;L ) across 206 mice. Cardiac gene expression was measured in 106 control (Pre-ISO) strains and 21 days after injection of isoproterenol (post-ISO) in 100 HMDP strains (a smaller number due to higher mortality of certain strains). Note that the significant correlation holds when considering separately pre-ISO (blue points, r ¼ 0:59, p ¼ 2 10 À11 ) and post-ISO (red points, r ¼ 0:42, p ¼ 1:5 10 À5 ) data. Lines show best fits of a linear model for pre-ISO (blue), post-ISO (red), and pre-and post-ISO combined (black). Expression data is taken from Santolini et al. (2018) and is averaged over all microarray probes for each gene. DOI: https://doi.org/10.7554/eLife.36717.016 also constraining the CaT. However, at least in the present study, constraining the AP waveform does not reproduce the correlation between Ca 2+ and K + currents observed in the HMDP.
Even though our calcium centric GES search did not include the conductance of the Na + current, it seems physiologically plausible that this conductance (and gap junction coupling) could also be regulated by Ca 2+ sensing to ensure that conduction velocity is adequate to generate a synchronous blood pressure waveform. In particular, the integrated CaT of the ventricles has to be reasonably synchronous to generate a normal blood pressure waveform, requiring the Na current density to be adequate for a normal conduction velocity through the tissue.
Even though the GES search was performed using target values of Ca 2+ sensors for a 4 Hz pacing frequency, different GES exhibit similar CaT amplitude versus pacing frequency curves consistent with experimental measurements (Figure 18 in Bondarenko et al., 2004) over a broad range of pacing frequencies from 0.5 to 4 Hz. Since different model parameters corresponding to different strains reproduce similar CaT amplitude-frequency curves, we do not expect the choice of pacing frequency to be critically important for calibrating model parameters that produce a normal electrophysiological phenotype. We attribute the robustness of those curves to the knock on effect of voltage on the L-type Ca 2+ current and SR Ca 2+ release via CICR. As a result of this effect, constraining the CaT indirectly constrains the relative magnitudes of depolarizing and repolarizing currents affecting the AP; that is, the same CaT amplitude can be obtained with combinations of Ca 2+ and K + currents that are both large or both small, thereby compensating each other, but not with combinations in which the Ca 2+ current is large and the sum of K + currents is small or vice-versa. It remains that the AP waveform and duration are only partially constrained by the CaT and are thus more variable than in a GES search that uses AP features such as duration, plateau voltage, etc., to constrain parameter sets (Sarkar and Sobie, 2009;Sarkar and Sobie, 2010).
While the additional constraint to keep the intracellular Na + concentration Na ½ i within a normal physiological range is not necessary to produce a physiological AP waveform, it constrains more tightly the conductance of the Na + -Ca 2+ exchanger current (compare I NaCa histograms of Figure 3D and Figure 3-figure supplement 1). This leads in turn to a tighter three-way compensation between the L-type Ca 2+ current and two dominant repolarizing K + currents in the mouse, which is less pronounced in a GES search in which Na ½ i is not constrained (Figure 3-figure supplement 2) and I Ca;L can be compensated by I NaCa in addition to those K + currents (four-way compensation). In the present study, we have leveraged the fact that different mice in the same HMDP inbred strain are isogenic to distinguish for the first time intra-heart cell-to-cell variability from inter-subject (interstrain in the HMPD context) variability. This has allowed us to use several hearts for each strain and perform current measurements in enough cells to statistically distinguish mean conductances of several currents in different strains. The results ( Figure 4A, Table 1 and Figure 4-figure supplement 1) clearly show that conductances differ between strains. Statistical testing shows that I Ca;L conductance measurements for different strains are very unlikely to belong to the same distribution (as indicated by the very small p-value) and mean conductances can vary by as much as two-hand-a-half fold between pairs of strains (e.g. C57BL/6J and BXA25/PgnJ), far in excess of the typical standard error of the mean. Furthermore, guided by the predictions of our computational GES search, we have also measured conductances of two dominant repolarizing currents in the mouse, I to;f and I Kur , to test for the existence of compensation between those currents and I Ca;L . The results ( Figure 4A) show that for eight out of the nine strains in which all three currents were measured, the currents accurately compensate each other as predicted by the GES search in which both the CaT and Na ½ i are constrained. Compensation is evidenced by the linear regression fit of I to;f +I Kur versus I Ca;L . One outlier strain (BXA12/PgnJ) deviates from this fit but still falls within the larger ensemble of computationally predicted GES without the Na ½ i constraint. Importantly, cells isolated from mouse strains with very different I Ca;L conductance have statistically indistinguishable contractile function ( Figure 4B). This suggests that compensation between I Ca;L and K + currents in different HMDP strains is present to maintain Ca 2+ homeostasis, as assumed in the computational GES search. As a whole, the results clearly support the hypothesis that Ca 2+ concentration plays a major role in feedback sensing of cellular activity and regulation of ion channel expression. While more strains would need to be studied to more accurately determine the role of the Na + concentration, measurements reported in Figure 4A suggest that it plays at least an auxiliary role in further constraining conductances beyond Ca 2+ sensing.
One limitation of the present GES search is that it only identifies possible combinations of electrophysiological parameters underlying a normal cardiac electrophysiological phenotype. However, it cannot by itself predict which GES among those found are actually represented, even approximately, in a genetically diverse population. Another limitation is that dominant K + currents and AP features are markedly different in mouse than human. While there is no analog of the HMDP for human, or other species (such as rabbit, dog, or pig) with AP features similar to human, it might be possible to extend the present study using renewable cardiomyocytes (CMs) derived from induced pluripotent stem cells (iPSC) as an alternative to the HMDP. However, iPSC-CMs and adult myocytes isolated from intact hearts exhibit quantitative differences in their responses to ionic current perturbations (Gong and Sobie, 2018). Therefore, it is unclear whether the variation of ionic conductances in iPSC-CMs of genetically different subjects would be representative of the variation of conductances in intact hearts of the same subjects, which is ultimately relevant for pharmacological treatment of cardiac arrhythmias. In addition, from the results of a recent study of AP variability in cardiomyocytes derived from different human subjects (Britton et al., 2017), it is unclear if inter-and intra-subject variability could be statistically distinguished in a large population.
Finally, the correlation between I Ca;L conductance and cardiac hypertrophic response of HMDP strains to sustained b-adrenergic stimulation ( Figure 6) also highlights the importance of considering the inherent variability of electrophysiological parameters in a genetically diverse population to interpret the variability of phenotypic response to pharmacological perturbations (Sarkar et al., 2012;Britton et al., 2013;Britton et al., 2017;Gong and Sobie, 2018) or stressors (Rau et al., 2017;Santolini et al., 2018). For example, pharmacological treatment with an anti-arrhythmic L-type calcium channel blocker, or pathologies such as hyperkalemia (elevated potassium level), would be expected to have different effects in different subjects. In the setting of the HMDP, an L-type calcium channel blocker or hyperkalemia would be expected to have a stronger effect on the calcium transient and action potential of mice strains that function under normal conditions with larger I Ca;L and potassium current conductances. Consistent with previous population level studies (Sarkar et al., 2012;Britton et al., 2013;Britton et al., 2017;Gong and Sobie, 2018;Rau et al., 2017;Santolini et al., 2018), taking into account this variability seems ultimately needed to develop personalized therapies for cardiac arrhythmias and heart failure. Table 1. Patch clamp measurements of I Ca;L , I to;f , and I Kur functional current density. Mean current density averaged over n cells isolated from multiple hearts for each strain is given together with the standard error.

Overview of the HMDP
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.

Electrophysiological and contraction measurements Cell isolation
Ventricular myocytes were enzymatically isolated from the hearts of adult female mice (8-to 12week-old) using a procedure previously developed and utilized to isolate rabbit cardiomyocytes by Yang et al. (2008). Briefly, hearts were removed from mice anesthetized with intravenous pentobarbital and perfused retrogradely at 37˚C in Langendorff fashion with nominally Ca 2+ -free Tyrode's buffer containing 1.2 mg/ml collagenase type II (catalog number 4176; Worthington) and 0.12 mg/ ml protease type XIV (catalog number P5147; Sigma) for 10 -17 min. After washing out the enzyme solution, the ventricles were cut from the atria and aorta and transferred to a separate glass dish containing Tyrode's solution. Cells were isolated by gentle mechanical dissociation, stored at room temperature, and used within 5 hr. This procedure typically yielded 30 -50% of rod-shaped and Ca 2+ -tolerant myocytes.

Patch clamping
Isolated ventricular myocytes were patch clamped in the whole cell ruptured patch configuration using borosilicate glass pipettes (1À3-megaohm tip resistance). Myocytes were superfused at 34À36˚C with Tyrode's solution modified accordingly. Currents were measured under voltage clamp conditions, using an Axopach 200B amplifier with a Digidata 1440A interface (Axon Instruments, Union City, CA). Data were acquired and analyzed using pClamp (Axon instruments) and Origin (Origin).

Potassium current measurements
For characterization of K + current properties, the pipette solution, designed to eliminate Ca 2+ currents, contained (in mM) 130 KCl, 5 MgATP, five phophocreatine di(tris), 5 HEPES, 5 NaCl, and 10 BAPTA (pH adjusted with HEPES to 7.1 -7.2).The superfusate, designed to eliminate Ca 2+ currents, contained (in mM) 136 NaCl, 5.4 KCl, 1 MgC12, 0.33 NaH2PO4, 10 glucose, 5 HEPES, and 0.2 CdCl2 (pH adjusted with Trizma base to 7.4). To access the activation of K + current components: I Kur , I to;f and I Kss , we adopted a voltage protocol, similar to the one reported by Zhou et al. (1998), in combination with the usage of 4-aminopyridine (4-AP) that has the following pharmacological properties: (1) I Kur is markedly blocked by 4-AP at submillimolar concentration (e.g. 0.1 mM); (2) a higher concentration (i.e., >1 mM) blocks I to;f effectively; and (3) I Kss is 4-AP resistant. Using this procedure, peak I Kur , I to;f and I Kss currents could be deduced from three current measurements without 4-AP and with 0.1 mM and 2 mM 4-AP after a voltage step from À50 to 0 mV.

Results of patch clamp measurements
Results of patch clamp measurements for all strains are summarized in Table 1 and Figure 4-figure supplement 1. Both the L-type Ca 2+ current and two K + currents (I to;f and I Kur ) were measured in nine strains and the L-type Ca 2+ current alone was measured in 16 strains.

Contraction analysis
The length of ventricular myocytes was measured and analyzed using the method described by Sdek et al. (2011). Briefly, myocytes were imaged during pacing using a high-speed charge-coupled device-based camera (128 Â 128 pixels; Cascade 128+; Photometrics) at 290 frames/s. The acquired video image data were then processed using Imaging Workbench software (version 6.0; INDEC Bio-Systems). Myocyte length was measured and analyzed using ImageJ software. Cell shortening was calculated from the ratio of peak systolic length to resting diastolic length averaged over 10 contractions evoked by the stimulus train. Results of cell shortening measurements for six strains paced at 4 Hz are summarized in Table 2.
Heart extraction and mass measurement for cardiac hypertrophic response 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).

Mathematical model of mouse ventricular myocytes
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 2+ spark recruitment rate to I Ca;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 (  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.

Equations for Ca 2+ cycling
We use the model for Ca 2+ cycling developed by Shiferaw et al. (2003) and subsequently implemented in Mahajan et al. (2008). The equations for Ca 2+ cycling are (2) where the SR leak flux and RyR release flux are given by

Intracellular Ca 2+ buffering
Similarly to Mahajan et al. (2008). All buffering parameters are experimentally based and summarized in Shannon et al. (2004). Buffering to SR, calmodulin, membrane, and sarcolemma binding sites are modeled using the instantaneous buffering approximation given by   Table 5 continued on next page Buffering to Troponin C is given by The SERCA uptake pump Similarly to Shiferaw et al. (2003). Na + dynamics Intracellular Na + dynamics are given by In order to reduce computation time, we have sped up the rate at which the system reaches steady-state by increasing d Na þ ½ i =dt by a factor of 100. This will make sodium converge to steadystate on a time-scale fast enough to perform the number of simulations necessary for this study. Once the cell reaches steady-state, d Na þ ½ i =dt is zero, so this modification will not affect the sodium dynamics at steady-state. By doing this, we can save up to 90% of calculation before the system reaches steady-state.

Ionic currents
The rate of change of the membrane voltage V is described by the equation where I stim is the external stimulus current driving the cell.
Calcium background leak (I Ca;b ) The sarcolemmal Ca 2+ ATPase (I PMCA ) The sarcolemmal Ca 2+ pump (I PMCA ) provides another mechanism, in addition to the exchanger (I NaCa ), for the extrusion of Ca 2+ ions out of the cell. This pump is not included in Mahajan et al. (2008). We added this current using the formula used by Bondarenko et al., 2004.
The Na + -Ca 2+ exchange flux (NaCa) The fast sodium current (I Na ) Similarly to the Shannon et al. (2004) model, Sodium background leak (I Na;b ) Inward rectifier K þ current (I K1 ) Similarly to the Bondarenko et al., 2004 model, where r ¼ 0:37 accounts for the presence of a persistent outward potassium current in patch clamp measurements of I to;f . We have increased the rates of the inactivation gate (a i and b i ) from the original formulation to match experimental measurements of I to;f inactivation rate under voltage clamp.

Computational reproduction of patch clamp experiments
In order to compare the GESs found by the computational search to the phenotype variability found in the HMDP, we iterate the model with voltage held constant, reproducing the experimental patch clamp procedure described in the Methods section of the main text. Each model is simulated for 1 s with V m held at À50 mV in order to reach steady-state. V m is then raised to 0 mV, and peak values I Ca;L , I to;f and I Kur are recorded.

Tissue scale modeling
Tissue scale modeling is performed using a 56 Â 56 array of myocytes, each with individual values of ionic conductances. Electrotonic coupling is simulated by introducing a diffusive into the V m evolution equation, The applied stimulus current occurs at a pacing rate of 4 Hz and is applied to each myocyte simultaneously. The diffusive term is applied isotropically, with diffusive co-efficient, D ¼ 1 cm 2 =s. In the discretized diffusion equation, We use a lattice size of Dx ¼ 225m, such that the 56 Â 56 lattice represents a 1.25 cm Â1.25 cm tissue.

GES search
In this study, we consider variation in six important ionic currents: L-type Ca current (I Ca;L ), the SR ATPase SERCA, Na + -Ca 2+ exchange (NaCa), ryanodine receptor (RyR), the transient outward K + current (I to;f ), and the ultra-rapidly-activating K + current (I Kur ). Those six currents were selected because they are the major currents influencing the CaT. The strength of these ionic currents is determined by their conductance g i . Any given set of parameters (p ¼ p 1 ; p 2 ; Á Á Á ; p n f g) corresponds to a different candidate myocyte model and produces a different phenotype, which we characterize by quantifiable measurements of its steady-state behaviour (sensors) that is steady state calcium transient amplitude, action potential duration and sarcoplasmic reticulum (SR) Ca 2+ concentration. When stimulating a model with a given period, these parameters (once the simulation has reached steady state) produce a phenotype which we can compare to the phenotype produced by the standard parameters of our model (p ref ).
We define a cost function E that quantifies how much each model's phenotype differs from our reference phenotype as where the S n p ð Þs are sensors characterizing the electrophysiological phenotype of the model's output. E p ref ð Þ is zero by definition. The three sensors used in this study are listed in Table 6. For any given set of conductances, a simulation is performed that outputs a value for each sensor. All values are calculated after the simulation is paced with a pacing cycle length PCL = 250 ms for 12.5 s when the system has reached steady state.
We define a good enough solution (GES) as a set of conductances with a phenotype such that the value of cost function E p ð Þ is less than a threshold ¼ 0:05. None of the cost function sensors S n are based on the membrane potential, and therefore a GES does not necessarily have an action potential shape close to the reference action potential shape. A GES is required to achieve steadystate. Therefore, parameters that produce parameter sets that do not reach a steady state during pacing at constant cycle length, such as those which exhibit calcium transient alternans, are not considered GESs. We additionally reject any set of parameters for which the output steady-state SR Ca 2+ load is above a threshold Ca 2þ Â Ã SR > 130 M Cyt , which we consider to be unphysiologically overloaded.