1. Computational and Systems Biology
  2. Physics of Living Systems
Download icon

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

  1. Colin M Rees
  2. Jun-Hai Yang
  3. Marc Santolini
  4. Aldons J Lusis
  5. James N Weiss
  6. Alain Karma  Is a corresponding author
  1. Northeastern University, United states
  2. Northeastern University, United States
  3. Cardiovascular Research Laboratory, David Geffen School of Medicine, University of California, United states
  4. David Geffen School of Medicine, University of California, United States
Research Article
  • Cited 0
  • Views 478
  • Annotations
Cite this article as: eLife 2018;7:e36717 doi: 10.7554/eLife.36717

Abstract

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.

https://doi.org/10.7554/eLife.36717.001

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 (Ca2+) 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+/Ca2+ exchanger, and Ca2+ handling proteins such as the ryanodine receptor (RyR) Ca2+ release channels, which open in response to Ca2+ entry into the cell via L-type Ca2+ channels, and the sarcoplasmic reticulum (SR) Ca2+ ATPase (SERCA), which uptakes Ca2+ back into the SR. Ca2+ release and uptake from the SR causes a transient rise in cytosolic Ca2+ 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 ‘species-averaged’) 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.

Schematic representation of the sarcolemmal currents and intracellular Ca2+ cycling proteins of the mouse ventricular myocyte model.
https://doi.org/10.7554/eLife.36717.002

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 Ca2+ 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 (Vm) 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 Vm and Ca2+ signals can further constrain model parameters (Sarkar and Sobie, 2009; Sarkar and Sobie, 2010). Groenendaal et al., 2015 constrained model parameters using Vm 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 Vm 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 (IKr 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 (Passini et al., 2016; 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 ‘Ca2+ 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 Ca2+, 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 Ca2+ 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 Ca2+/calcineurin/NFAT pathway regulates L-type Ca2+ 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 Ca2+ current (ICa,L) and K+ currents (Ito,f and IKur) 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 cell-to-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 (ICa,L, Ito,f and IKur) 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 (ICa,L current is large or small when the sum of Ito,f and IKur currents is large or small, respectively) predicted by our computational Ca2+ 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).

Results

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 Δ[Ca]i between the diastolic and peak value of the cytosolic Ca2+ concentration [Ca]i, and the time-averaged value of [Ca]i over one pacing period, denoted by [Ca]i. The CaT amplitude Δ[Ca]i is a major determinant of the contractile force while [Ca]i provides an average measure of the cytosolic Ca2+ concentration in the cell. Both quantities will be used as Ca2+ 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 Ca2+ handling proteins that include the ryanodine receptor (RyR) Ca2+ release channels and the sarcoplasmic reticulum (SR) Ca2+ ATPase SERCA, which pumps Ca2+ 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 Δ[Ca]i and [Ca]i and the intracellular sodium concentration [Na]i no longer vary from beat to beat.

Effects of individual conductances on the Ca2+ transient (CaT).

(A) CaT amplitude defined as the difference Δ[Ca]i between the peak and diastolic values of the cytosolic Ca2+ concentration [Ca]i versus G/Gref where G is the individual conductance value and Gref some fixed reference value. (B) Time-averaged [Ca]i over one pacing period ([Ca]i) versus G/Gref. Illustration of the effect of varying ICa,L conductance (C) and Ito,f conductance (D) on AP and CaT profiles, where 50%, 100%, and 150% correspond to Gref=0.5, 1.0, and 1.5, respectively. (E) Effect of varying RyR conductance on SR Ca2+concentration [Ca]SR and CaT. Different time windows are plotted for the CaT and SR load (0 to 150 ms) and AP waveforms (0 to 50 ms) in (C-–E).

https://doi.org/10.7554/eLife.36717.003

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/Gref where Gref is a reference value producing a normal CaT. Increasing the conductance of the inward L-type Ca2+ current ICa,L is seen to strongly increase both Δ[Ca]i and [Ca]i but has a weak effect on the AP waveform (Figure 2C). The effect on the CaT stem from the fact that ICa,L is the main trigger of Ca2+-induced Ca2+ release(CICR), which transfers a large amount of Ca2+ from the SR to the cytosol. The weak effect on APD is due to the fact that the increase of CaT amplitude causes ICa,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 Ito,f (fast inactivating component of the transient outward current) and IKur causes both Δ[Ca]i and [Ca]i to decrease. Increasing either of those K+ currents speeds up repolarization (as illustrated for Ito,f in Figure 2D) and hence inactivation of ICa,L, thereby reducing the magnitude of CICR. Increasing the conductance of the sodium-calcium exchanger current INaCa also causes both Δ[Ca]i and [Ca]i to decrease by enhancing the forward mode of this current that extrudes Ca2+ 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 ICa,L conductance on CaT amplitude in the present mouse model. Changing RyR expression from its reference value is seen to leave Δ[Ca]i and [Ca]i almost unchanged, even though it strongly affects the SR Ca2+ 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 Ca2+ 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 Δ[Ca]i and [Ca]i. Increasing SERCA magnitude increases SR Ca2+ load, thereby increasing the amount of SR Ca2+ release and CaT amplitude, but at the same time depletes Ca2+ 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 Δ[Ca]i, time averaged cytosolic Ca2+ concentration [Ca]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 Ca2+ sensors, instead of Δ[Ca]i (the difference between the peak and diastolic [Ca]i values) and time averaged [Ca]i, gives nearly identical results. So our Ca2+ 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 Δ[Ca]i, [Ca]i, and [Na]i corresponding to the reference set of parameters (Gref values) of the ventricular mycoyte model. The search was conducted by defining a cost function

(1) E(p)=n=1N (Sn(p)SnSn)2ϵ,

which is an aggregate measure of the deviation of output sensors Sn(p) from their desired target values Sn. Here, N=3 with S1=Δ[Ca]i, S2=[Ca]i, and S3=[Na]i, and ϵ is a small tolerance that we choose to be 5%. E is a function of model parameters p=(p1, p2,) chosen to consist of the conductances of ICa,L, Ito,f, IKur, and INaCa as well as RyR and SERCA expression levels. Effects of individual changes of those parameters on CaT properties measured by S1 and S2 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 (p1, p2,) a value comprised between 0% and 300% of its reference value Gref. 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.

Figure 3 with 3 supplements see all
Computationally determined good enough solutions (GES) with calcium sensing.

(A) Examples of GES representing combinations of 6 conductances that produce a normal CaT and intracellular Na+ concentration. Each color represents a different GES and the corresponding AP and CaT profiles are shown in B) and C), respectively. (D) Histograms of individual normalized conductances G/Gref for a collection of 7263 GES showing that some conductances are highly variable while others are highly constrained. (E) Three-dimensional (3D) plot revealing a three-way compensation between conductances of ICa,L, Ito,f, and IKur. 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 ICa,L versus the sum of Ito,f and IKur. 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.

https://doi.org/10.7554/eLife.36717.004

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 INaCa, which turns out to be constrained by the intracellular sodium concentration sensor (S3=[Na]i). This is revealed by a GES search with only Ca2+ sensing (S1 and S2) that yields a broader histogram for the INaCa 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ 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 ICa,LIto,f and IKur. 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 ICa,L, Ito,f, and IKur 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 Ca2+ protein expression levels). However, SERCA expression and INaCa conductance are constrained by Ca2+ and Na+ sensing, respectively, and RyR expression has a negligible effect on both Ca2+ 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 INaCa 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 ICa,L, Ito,f, and IKur conductances by plotting the sum of the peak currents of Ito,f and IKur versus the peak current of ICa,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 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 ICa,L, Ito,f, and IKur 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 Ca2+ 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 Ca2+ 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 Ito,f and IKur versus the peak current of ICa,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 Ca2+ 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 Ito,f and IKur follows a linear correlation with ICa,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 ICa,L measurements. The result shows that ICa,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 ICa,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 Ca2+ and K+ currents were measured for the nine strains reported in Figure 4A, the Ca2+ 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 ICa,L conductances.

Figure 4 with 1 supplement see all
Good enough solutions in the Hybrid Mouse Diversity Panel (HMDP).

(A) Central result of this paper showing quantitative agreement between theoretically predicted and experimentally measured compensation of inward Ca2+ and outward K+ currents. Equivalent plot of Figure 3F showing the sum of Ito,f and IKur versus ICa,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 Ca2+ and K+ currents. The three strains selected for the organ scale study (C57BL/6J, CXB1/ByJ, and BXA25/PgnJ) with low, medium, and high ICa,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.

https://doi.org/10.7554/eLife.36717.008
Table 1
Patch clamp measurements of ICa,L, Ito,f, and IKur functional current density.

Mean current density averaged over n cells isolated from multiple hearts for each strain is given together with the standard error.

https://doi.org/10.7554/eLife.36717.010
StrainICa,L (pA/pF)nIKur (pA/pF)nIto,f (pA/pF)nIKss (pA/pF)n
A/J13.71947 ± 1.230851911.61804 ± 2.790891313.03735 ± 1.64818145.23445 ± 0.4837713
BALB/cByJ10.72278 ± 1.3951911.496 ± 2.03274107.46938 ± 0.8261598.197 ± 0.6243810
BTBR T+tf/J14.75667 ± 1.14159616.42364 ± 2.78295119.271 ± 1.43985106.093 ± 0.395539
BXA12/PgnJ11.01333 ± 0.98995915.66429 ± 1.70548811.89875 ± 2.1802788.34556 ± 1.300659
BXA25/PgnJ17.57 ± 4.81376515.8049 ± 1.74509711.62704 ± 1.5961377.11986 ± 1.315546
BXH6/TyJ7.32625 ± 0.59327165.954 ± 0.43731108.73172 ± 0.67617117.98545 ± 0.8754511
C57BL/6J7.3925 ± 0.93181107.82727 ± 1.45134117.594 ± 0.92097109.16273 ± 1.6214811
CXB1/ByJ11.93909 ± 0.810221111.69556 ± 1.4209598.73883 ± 1.11087106.089 ± 0.572410
CXB11/HiAJ7.63625 ± 0.8934488.90316 ± 2.4975568.42537 ± 1.2463576.8481 ± 1.36387
AXB8/PgnJ7.355 ± 1.096126---
BXA14/PgnJ11.28091 ± 1.0079615---
BXA4/PgnJ11.77615 ± 1.1947613---
BXD34/TyJ11.12556 ± 1.137419---
CBA/J7.956 ± 0.8126910---
CXB7/ByJ9.05386 ± 0.608716---
SJL/J11.2745 ± 0.9970820---

Unlike ion channel conductances, CaT properties were assumed not to vary in the computationally-enabled GES search, which rests on the hypothesis that Ca2+ 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
Organ scale compensation.

(A) Mean ICa,L conductance in three different HMDP strains where thick and thin bars denote standard error and standard deviation, respectively. (B) Sets of conductances generated to be representative of individual cells within ventricular tissue of the three strains by assigning normally distributed random values to the ICa,L, Ito,f and IKur conductances using experimentally determined means and standard deviations. The blue, green, and red points correspond to the three HMDP strains with low (C57BL/6J), medium (CXB1/ByJ), and high (BXA25/PgnJ) ICa,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 Ca2+ transient (CaT) amplitude (ΔCa) 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 Ca2+ and K+ currents remains operative at a tissue scale. (G) Distribution of CaT amplitudes obtained by varying only ICa,L conductance and with Ito,f and IKur conductances fixed to their reference values. Lack of compensation between Ca2+ and K+ currents in this case yields different mean CaT amplitude.

https://doi.org/10.7554/eLife.36717.011

Compensation at the organ scale

Current measurements discussed in the previous section show that mean conductances of Ca2+ 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 ICa,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 Ca2+ 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 ICa,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 ICa,L, Ito,f, and IKur were assumed to vary randomly from cell to cell, and no constraint was imposed on the ratio of Ito,f+IKur to ICa,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 Ito,f+IKur to ICa,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 Δ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 Ito,f+IKur to ICa,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 ICa,L conductance (CXB1/ByJ), the different ICa,L conductances are not compensated by different Ito,f and IKur 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 ICa,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 (mpre) and 3 weeks after (mpost) implantation of a pump continuously delivering isoproterenol (Table 3). Figure 6 reveals the existence of a statistically very significant correlation between baseline ICa,L conductance and the hypertrophic response (mpost/mpre). Although many factors contribute to the hypertrophic response in the HMDP (Rau et al., 2017; Santolini et al., 2018), intracellular Ca2+ overload activating the Ca2+-calcineurin-NFAT signaling pathway has been shown to play a major role (Bers, 2008). Since ICa,L is the major pathway of Ca2+ entry into the cytoplasm, it is intriguing to speculate that strains with a larger baseline ICa,L conductance under baseline conditions have a more robust increase in ICa,L that is not adequately compensated by repolarizing K+ currents, making those strains more susceptible to Ca2+ overload when ICa,L is enhanced during sustained β-adrenergic stimulation by isoproterenol. Hypothetically, this may result in a stronger cardiac hypertrophic response. To make this case convincingly, however, would require demonstrating that Ca2+ overload is chronically worsened in strains with a high baseline ICa,L conductance and ruling out other strain-dependent hypertrophy-promoting pathways that are not Ca2+-dependent, which is beyond the scope of the present work.

Correlation between L-type Ca2+ current conductance and cardiac hypertrophic response to a stressor for different HMDP strains.

The Pearson correlation is r = 0.86 (p=3e-4).

https://doi.org/10.7554/eLife.36717.014
Table 3
Heart mass before and 3 weeks after Isoproterenol (ISO) injection.
https://doi.org/10.7554/eLife.36717.015
StrainHeart mass pre-ISO, mpre (g)Heart mass post-ISO, mpost (g)
A/J0.0886666670.133
AXB8/PgnJ0.0870.0992
BALB/cByJ0.101050.1389
BTBR T+tf/J0.141620.223
BXA-12/PgnJ0.064NA
BXA-14/PgnJ0.09750.1248
BXA-4/PgnJ0.10310.14675
BXD-34/TyJ0.1215NA
BXH-6/TyJ0.08450.0878
C57BL/6J0.0967166670.1222
CBA/J0.0953333330.13
CXB-11/HiAJ0.11350.1255
CXB-7/ByJ0.1090.1475
SJL/J0.0870.123333333

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 ICa,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 ICa,L, Ito,f, and IKur, respectively. However, a statistically very significant correlation was found between expression levels of Cacna1c and Kcnip2 that encodes the KChIP2 accessory β 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 Ito,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 ICa,L and Ito,f+IKur functional current densities (Figure 4).

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.11013) between the expression level of Kcnip2, encoding the KChIP2 accessory β subunits that interact with Kv4.2 channels (Ito,f) and of Cacna1c, a gene encoding the α1C subunit of the Cav1.2 L-type calcium channels (ICa,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.59p=21011) and post-ISO (red points, r=0.42p=1.5105) 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.

https://doi.org/10.7554/eLife.36717.016

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 Ca2+ 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 post-translational 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 Ca2+ (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 Ca2+ 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 (Xiao et al., 2008) has provided evidence of a compensatory increase of IKs following exposure of canine cardiomyocytes to a pharmacological IKr 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 IKs 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 Ca2+ sensing hypothesis and suggests its generality beyond mice.

A remarkable and nontrivial finding of the present computational study is that Ca2+ 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 Ca2+ 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 Ca2+ and outward K+ currents observed in the HMDP (Figure 4A) was no longer preserved, since other inward and outward currents including INaCa 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), Ca2+ 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 Ca2+ 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 also constraining the CaT. However, at least in the present study, constraining the AP waveform does not reproduce the correlation between Ca2+ 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 Ca2+ 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 Ca2+ 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 Ca2+ current and SR Ca2+ 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 Ca2+ and K+ currents that are both large or both small, thereby compensating each other, but not with combinations in which the Ca2+ 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+-Ca2+ exchanger current (compare INaCa histograms of Figure 3D and Figure 3—figure supplement 1). This leads in turn to a tighter three-way compensation between the L-type Ca2+ 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 ICa,L can be compensated by INaCa 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 (inter-strain 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 ICa,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, Ito,f and IKur, to test for the existence of compensation between those currents and ICa,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 Ito,f+IKur versus ICa,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 ICa,L conductance have statistically indistinguishable contractile function (Figure 4B). This suggests that compensation between ICa,L and K+ currents in different HMDP strains is present to maintain Ca2+ homeostasis, as assumed in the computational GES search. As a whole, the results clearly support the hypothesis that Ca2+ 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 Ca2+ 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 ICa,L conductance and cardiac hypertrophic response of HMDP strains to sustained β-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 ICa,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.

Materials and methods

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 12-week-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 Ca2+-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 Ca2+-tolerant myocytes.

Patch clamping

Isolated ventricular myocytes were patch clamped in the whole cell ruptured patch configuration using borosilicate glass pipettes (13-megaohm tip resistance). Myocytes were superfused at 3436°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).

L-type calcium current measurements

For characterization of ICa properties, the pipette solution, designed to eliminate K+ and Cl currents, contained (in mM) 100 CsMeS, 30 CsCl, 5 MgATP, five phophocreatine di(tris), 5 N-2-hydroxyethylpi-perazine-N’−2-ethanesulfonic acid (HEPES), 5 NaCl and 0.1 1,2-Bis(2-Aminophenoxy)ethane-N,N,N’,N’-tetraacetic acid (BAPTA) (pH adjusted with HEPES to 7.1 – 7.2). The superfusate, designed to eliminate K+ currents, contained (in mM), 136 NaCl, 5.4 CsCl, 1 MgC12, 0.33 NaH2PO4, 10 glucose, 5 HEPES, and 1.2 CaC12 (pH adjusted with Trizma base to 7.4). Voltage-clamp protocols to assess activation were as previously described by Delbridge et al. (1997). In all experiments, peak ICa current was recorded after a voltage step from −50 to 0 mV.

Potassium current measurements

For characterization of K+ current properties, the pipette solution, designed to eliminate Ca2+ 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 Ca2+ 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: IKur, Ito,f and IKss, 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) IKur is markedly blocked by 4-AP at submillimolar concentration (e.g. 0.1 mM); (2) a higher concentration (i.e., >1 mM) blocks Ito,f effectively; and (3) IKss is 4-AP resistant. Using this procedure, peak IKur, Ito,f and IKss 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 Ca2+ current and two K+ currents (Ito,f and IKur) were measured in nine strains and the L-type Ca2+ 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 BioSystems). 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.

Table 2
Cell Shortening at 4 Hz pacing.
https://doi.org/10.7554/eLife.36717.017
StrainΔL/L (4 Hz)n
BXA12/PgnJ0.0835 ± 0.01845
BXA14/PgnJ0.1182 ± 0.010814
BTBR T+tf/J0.0978 ± 0.00957
BALB/cByJ0.105 ± 0.01585
C57BL/6J0.0989 ± 0.0166
BXA25/PgnJ0.0838 ± 0.01125

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 Ca2+ spark recruitment rate to ICa,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 ICa,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 Ca2+ transient alternans. Even though the Shiferaw et al. (2003) model of calcium cycling model was developed for rabbit mycoytes, it can in principle describe Ca2+ 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 ICa,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 INa and ICa,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 Ca2+ cycling. We verified that the Shiferaw et al. (2003) model of Ca2+ cycling integrated with the Shannon et al. (2004) model of ICa,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.

Table 4
Reference values of ionic current parameters varied in the GES search.
https://doi.org/10.7554/eLife.36717.018
ParameterDefinitionReference valueReference source

gCa

Ca2+ current flux333.32 mmol/(Cm C)Measured

vup

Peak uptake rate1.17 μM/msChosen*

gNaCa

Peak NaCa rate36.6 μM/sBondarenko et al., 2004

gRyR

Release current
strength
12.9 sparks cm2/mAMahajan et al. (2008)

gto,f

Ito,f peak
conductance
0.16 A/FMeasured

gKur

IKur peak
conductance
0.144 A/FMeasured
  1. *vup was chosen such that the reference Ca2+ transient amplitude was normal.

Table 5
Mouse ventricular myocyte model parameters.
https://doi.org/10.7554/eLife.36717.019
ParameterDefinitionValue
Physical constants and ionic concentrations

Cm

Cell capacitance3.1 ×104μF

vi

Cell volume2.58 ×105μl

vs

Submembrane volume0.02 vi
FFaraday Constant96.485 C/mmol
RUniversal gas constant8.314 J mol1 K1
TTemperature298 K

[Na]o

External Na+
concentration
140 mM

[K]i

Internal K+
concentration
143.5 mM

[K]o

External K+
concentration
5.4 mM
[Ca2+]oExternal Ca2+
concentration
1.8 mM
Cytosolic buffering parameters

BT

Troponin C
concentration
70 μmol/l cyt

konT

on rate for Troponin C
binding
0.0327 (μM ms)1

koffT

off rate for Troponin C
binding
0.0196 (ms)1

BSR

SR binding site
concentration
47 μmol/l cyt

KSR

SR binding site disassociation constant0.6 μM

BCd

Calmodulin binding site concentration24 μmol/l cyt

KCd

Calmodulin binding site disassociation constant7 μM

Bmem

Membrane binding site concentration15 μmol/l cyt

Kmem

Membrane binding site disassociation constant0.3 μM

Bsar

Sarcolemma binding
site concentration
42 μmol/l cyt

Ksar

Sarcolemma binding
site disassociation constant
13 μM
SR release parameters

τr

Spark lifetime10 ms

τa

NSR-JSR diffusion time20 ms
 uRelease slope4 ms1

csr

Release slope threshold90 μM / l cytosol

τd

cp - cs diffusion time0.50 ms*

τs

cs - ci diffusion time0.75 ms
Exchanger, uptake, and SR leak parameters

cup

Uptake threshold0.5 μM

ksat

NaCa saturation
threshold
0.1

ξ

NaCa energy barrier
position
0.35

Km,Nai

Ion mobility constant21 mM

Km,Nao

Ion mobility constant87.5 mM

Km,Cao

Ion mobility constant1380 μM

gl

Leak current
conductance
1.74 × 105ms1
Ionic current parameters

gNa

Na+ current
conductance
13 mS/μF

gNa,b

Na+ background
current conductance
0.0026 mS/μF

gCa,b

Ca2+ background
current conductance
0.000367 mS/μ F

g¯Ca

Strength of local LCC
calcium flux
9000 mM/(cm C)

gK1

IK1 conductance0.2938 mS/μF

gNaK

INaK conductance1.716 mS/μF

gKss

IKss conductance0.025 mS/μF

gto,s

Ito,s conductance0 mS/μF

JPMCA,max

Maximal JPMCA fluxone pA/pF

KPMCA

Saturation constant for Ca2+ current0.5 μM

PCa

Constant0.00054 cm/s

Po,max

Constant0.083
  1. *We have reduced this value from the original value of Mahajan et al. (2008) so that the Ca2+ transient increases when SERCA uptake rate is increased.

Equations for Ca2+ cycling

We use the model for Ca2+ cycling developed by Shiferaw et al. (2003) and subsequently implemented in Mahajan et al. (2008). The equations for Ca2+ cycling are

(2) dcsdt=βs[vivs(JrelJd+JCa+JNaCa)Jtrpns],
(3) dcidt=βi[JdJup+JleakJtrpniJPMCA],
(4) dcjdt=Jrel+JupJleak,
(5) dcjdt=cjcjτa,
(6) Jd=csciτs,
(7) dcpdt=J~SR+J~Ca(cpcs)τd,

where the SR leak flux and RyR release flux are given by

(8) Jleak=gl(12.4 cjci),
(9) dJreldt=Ns(t)cjQ(cj)csrJrelT,
(10) T=τr1τr(dcjdt/cj)
(11) Q(cj)={ 0 0<cj<50,cj5050cjcsr,ucj+(1u)csr50cj>csr,
(12) Ns=gRyR(V)PoiCa,
(13) gRyR(V)=gRyRe0.05(V+30)1+e0.05(V+30),
(14) J~SR=gSR(V)Q(cj)PoiCa,
(15) gSR(V)=50 gRyR(V)

Intracellular Ca2+ 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

(16) βi=(1+BSRKSR(ci+KSR)2+BcdKcd(ci+Kcd)2+BmemKmem(ci+Kmem)2+BsarKsar(ci+Ksar)2)1,
(17) βs=(1+BSRKSR(cs+KSR)2+BcdKcd(cs+Kcd)2+BmemKmem(cs+Kmem)2+BsarKsar(cs+Ksar)2)1

Buffering to Troponin C is given by

(18) d[CaT]idt=Jtrpni=konTci(BT[CaT]i)koffT[CaT]i,
(19) d[CaT]sdt=Jtrpns=konTcs(BT[CaT]s)koffT[CaT]s

The SERCA uptake pump

Similarly to Shiferaw et al. (2003).

(20) Jup=vupci2ci2+cup2

Na+dynamics

Intracellular Na+ dynamics are given by

(21) d[Na+]idt=100CmFvi(INa+3INaCa+3INaK+INa,b)

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 steady-state 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

(22) dVdt=(Istim+ICa,L+IPMCA+INaCa+INa+Ito,f+Ito,s+IKur+IKss+IK1+INaK+ICa,b+INa,b)
(23) ICa,L=2FviCmJCa
(24) INaCa=FviCmJNaCa
(25) IPMCA=2FviCmJPMCA

where Istim is the external stimulus current driving the cell.

The L-type Ca current (ICa,L)

Similarly to the Shannon et al. (2004) model,

(26) JCa=gCaPoiCa
(27) J~Ca=g¯CaPoiCa
(28) iCa=4PCaVF2RTcse2VF/RT0.341[Ca2+]oe2VF/RT1
(29) Po=Po,max×d×f×fCa
(30) dfCadt=0.12(1fCa)1.40251+(30/cp)4fCa
(31) d=11+e(V+4.6)/6.3
(32) τd=d1e(V+4.6)/6.30.035(V+4.6)
(33) f=111+e(V+22.8)/6.1
(34) τf=10.020.007e(0.0337(V+10.5))2
(35) dddt=ddτd
(36) dfdt=ffτf

Calcium background leak (ICa,b)

(37) ICa,b=gCa,b(VECa)
(38) ECa=RT2Flog([Ca2+]oci)

The sarcolemmal Ca2+ ATPase (IPMCA)

The sarcolemmal Ca2+ pump (IPMCA) provides another mechanism, in addition to the exchanger (INaCa), for the extrusion of Ca2+ 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.

(39) JPMCA=JPMCA,maxci2KPMCA2+ci2

The Na+-Ca2+ exchange flux (NaCa)

Similarly to the Bondarenko et al., 2004 model,

(40) JNaCa=gNaCaeξVF/RT[Na+]i3[Ca2+]oe(ξ1)VF/RT[Na+]o3ci(1+ksate(ξ1)VF/RT)(Km,Nao3+[Na+]o3)(Km,Cao+[Ca2+]o)

The fast sodium current (INa)

Similarly to the Shannon et al. (2004) model,

(41) INa=gNam3hj(VENa)
(42) dhdt=αh(1h)βhh
(43) djdt=αj(1j)βjj
(44) dmdt=αm(1m)βmm
(45) αm=0.32V+47.131e0.1(V+47.13)
(46) βm=0.08eV/11
For V40mV,
(47) αh=0
(48) αj=0
(49) βh=10.13(1+e(V+10.66)/11.1)
(50) βj=0.3e2.535×107V1+e0.1(V+32)
For V40mV,
(51) αh=0.135e(V+80)/6.8
(52) βh=3.56e0.079V+3.1×105e0.35V
(53) αj=(1.2714×105e0.2444V3.474×105e0.04391V)×(V+37.78)1+e0.311(V+79.23)
(54) βj=0.1212e0.01052V1+e0.1378(V+40.14)

Sodium background leak (INa,b)

(55) INa,b=gNa,b(VENa)
(56) ENa=RTFlog([Na+]o[Na+]i)

Inward rectifier K+ current (IK1)

Similarly to the Bondarenko et al., 2004 model,

(57) IK1=gK1[K+]o[K+]o+0.21(VEK1+e0.0896(VEK))
(58) EK=RTFlog([K+]o[K+]i)

The fast component of the transient outward K+ current (Ito,f)

This current is modified from the formulation of Bondarenko et al., 2004 as:

(59) Ito, f=gto, fato, f3ito, f(VEK)
(60) dato, fdt=aato, fτa
(61) dito, fdt=iito, fτi
(62) αa=0.18264e0.03577(V+45)
(63) βa=0.3956e0.06237(V+45)
(64) αi=0.00152e(V+13.5)/7.00.067083e(V+33.5)/7.0+1
(65) βi=0.0095e(V+33.5)/7.00.051335e(V+33.5)/7.0+1
(66) τa=1αa+βa
(67) τi=1αi+βi
(68) a=αaαa+βa
(69) i=(1r)αiαi+βi+r

where r=0.37 accounts for the presence of a persistent outward potassium current in patch clamp measurements of Ito,f. We have increased the rates of the inactivation gate (αi and βi) from the original formulation to match experimental measurements of Ito,f inactivation rate under voltage clamp.

The slow component of the transient outward K+ current (Ito,s)

Similarly to the Bondarenko et al., 2004 model,

(70) Ito, s=gto, sato, sito, s(VEK)
(71) dato, sdt=assato, sτta, s
(72) dito, sdt=issito, sτti, s
(73) ass=11+e(V+22.5)/7.7
(74) iss=11+e(V+45.2)/5.7
(75) τta, s=0.493e0.0629V+2.058
(76) τti, s=270.0+10501+e(V+45.2)/5.7

The ultra-rapidly activating component of the delayed rectifier K+ current (IKur/IK,slow)

Similarly to the Bondarenko et al., 2004 model,

(77) IKur=gKurauriur(VEK)
(78) daurdt=assaurτaur
(79) diurdt=issiurτiur
(80) τaur=0.493e0.0629V+2.058
(81) τiur=12001701+e(V+45.2)/5.7

We have reduced the timescale of inactivation (τiur) from the original formulation to match experimental measurements of IKur inactivation rate under voltage clamp.

The non-inactivating steady-state K+ current (IKss)

Similarly to the Bondarenko et al., 2004 model,

(82) IKss=gKssaKss(VEK)
(83) daKssdt=assaKssτKss
(84) τKss=39.3e0.0862V+13.17

The Na+-K+ pump current (INaK)

Similarly to the Bondarenko et al., 2004 model,

(85) INaK=gNaKfNaK11+(Km,Nai/[Na+]i)3/2[K+]o[K+]o+Km,Ko
(86) fNaK=11+0.1245e0.1VF/RT+0.01548767σeVF/RT
(87) σ=17(e[Na+]o/673001)

Effect of cs-ci diffusion rate on how the Ca2+ transient depends on SERCA

The effect of modifications of the SERCA pump on the steady-state Ca2+ transient is shown in Figure 1. While increasing SERCA peak uptake current has the effect of sequestering Ca2+ back into the SR which would reduce the Ca2+ transient, the dominant effect is to increase the Ca2+ transient due to higher SR Ca2+ load at steady-state. This is consistent with experiments showing restoration of Ca2+ transient amplitude when SERCA is up regulated (del Monte et al., 2002del Monte et al., 1999).

Initial simulations using the original Ca2+ cycling parameters of Mahajan et al. (2008) showed that the SR Ca2+ load decreased as the uptake rate was increased since more Ca2+ was extruded from the sub membrane region of the myocyte via NCX before it had sufficient time to diffuse from the sub membrane space (with local Ca2+ concentration denoted by cs) into the cytosol compartment (with local Ca2+ concentration denoted by ci) to be re-uptaken into the SR. In order to rectify this, we increased the diffusion rate between the sub membrane and cytosolic compartments. We found that when this rate is faster (smaller τs), the SR Ca2+ load increases with increasing SERCA uptake rate as experimentally observed. For this reason, we use a value of τs = 0.75 ms for all simulations in this study.

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 Vm held at −50 mV in order to reach steady-state. Vm is then raised to 0 mV, and peak values ICa,L, Ito,f and IKur 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 Vm evolution equation,

(88) Vmt=1Cm(Istim+Iion)+D2Vm

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=1cm2/s. In the discretized diffusion equation,

(89) Vi, jt+1=Vi, jt+dt[1Cm(Istim+Iion)+DΔx2(Vi±1,j±1t4Vi, jt)],

We use a lattice size of Δx=225μm, 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 (ICa,L), the SR ATPase SERCA, Na+-Ca2+ exchange (NaCa), ryanodine receptor (RyR), the transient outward K+ current (Ito,f), and the ultra-rapidly-activating K+ current (IKur). 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 gi. Any given set of parameters (p={p1, p2,,pn}) 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) Ca2+ 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 (pref).

We define a cost function E that quantifies how much each model’s phenotype differs from our reference phenotype as

(90) E(p)2=n=1N(Sn(p)Sn(pref)Sn(pref))2,

where the Sn(p)s are sensors characterizing the electrophysiological phenotype of the model’s output. E(pref) is zero by definition. The three sensors used in this study are listed in Table 6.

Table 6
Simulation outputs corresponding to physiological sensors
https://doi.org/10.7554/eLife.36717.020
AbbreviationDescriptionReference value

[Ca]i

Average cytostolic Ca2+ over one beat0.24 μM

Δ[Ca]i

Ca2+ transient amplitude0.5 μM

[Na]i

Diastolic Na+14 mM

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 Sn 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 steady-state. 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 Ca2+ load is above a threshold [Ca2+]SR> 130 μMCyt, which we consider to be unphysiologically overloaded.

Completing an exhaustive search of the parameter space becomes increasingly computationally intensive as the number of parameters grows. An exhaustive search of a M dimensional space, considering K possible values for each conductance requires KM evaluations of the cost function. Additionally, as the number of sensors (N) that we use to calculate the cost function increases, the fraction of models tested that are good enough solutions (Φ) will decrease. We calculated E(p) for 107 random parameter sets, p, and found only 11 which were GESs (E(p)<0.05). It is therefore apparent that an exhaustive search is not efficient for finding GES in high dimensional parameter space and for this reason we use a minimization scheme to find GESs. We start by randomly assigning values to each parameter such that 0<pi<3pi,ref for each i in p, and then minimize them with respect to the cost function, E(p), using the Nelder-Mead simplex algorithm (Nelder and Mead, 1965) (also known as the Amoeba algorithm) until E(p)<ϵ. Running this procudure 10,000 times yeilds 7263 GESs, with the remaining trials being rejected either because the minimization algorithm does not converge, the SR load constraint is not satisfied, or the system is found to not be in steady-state (determined by comparing Δ[Ca]i of the 50th and 51 st beat).

The Nelder-Mead simplex algorithm

The Nelder-Mead algorithm (Nelder and Mead, 1965) maintains a non-degenerate simplex at each iteration, a geometric figure in n dimensions of nonzero volume that is the convex hull of n+1 vertices, x0,x1,...,xn, and their respective function values. Suppose we start from the vector x0, the simplex can be initialized as xi=x0+δei, where ei is a unit vector, and where δ is our guess of the problem’s characteristic length scale. In each iteration, new points are computed, along with their function values, to form a new simplex. The algorithm terminates when the function values at the vertices of the simplex satisfy a predetermined condition. One iteration of the Amoeba algorithm consists of the following steps (the standard values for the coefficients are: α=1,β=2,γ=0.5,σ=0.5):

  1. Order: order and re-label the n+1 vertices as x0,x1,...,xn, such that F(x0)F(x1)...F(xn). Since we want to minimize F(x0), we refer to x0 as the best point, to xn as the worst point, and to xn1 as the next worst point. Let xc refer to the centroid of the n points in the vertex. xc=i=0n1 xi/n.

  2. Reflect: compute the reflected point, xr=xc+α(xcxn). Evaluate F(xr). If F(x0)F(xr)<F(xn1), then obtain a new simplex by replacing the worst point xn with the reflected point xr and go to step 1.

  3. Expand: if F(xr)<F(x0), compute the expanded point, xe=xc+β(xrxc). If F(xe)<F(xr), then obtain a new simplex by replacing the worst point xn with the expanded point xe and go to step 1; otherwise then obtain a new simplex by replacing the worst point xn with the reflected point xr and go to step 1.

  4. Contract: At this step, where it is certain that F(xr)>F(xn1), compute the contracted point xcon=xc+γ(xnxc). If F(xcon)F(xn), obtain a new simplex by replacing the worst point xn with the expansion point xcon and go to step 1.

  5. Shrink: replace all vertices except the best x0 with xi=x0+σ(xix0) and go to step 1.

Two sensor search only constraining the Ca2+ transient

The results of the GES search described in the main text were reproduced using only two sensors, constraining [Ca]i and Δ[Ca]i but not constraining [Na]i. This results in a broader histogram for the INaCa conductance (Figure 3—figure supplement 1 compared to Figure 3).

Code availability

The codes used in this work are available at: https://github.com/circs/GES (Rees, 2018; copy archived at https://github.com/elifesciences-publications/GES).

References

  1. 1
  2. 2
  3. 3
    Excitation-Contraction Coupling and Cardiac Contractile Force, 237
    1. D Bers
    (2001)
    Springer Science & Business Media.
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
    Cardiac myocyte volume, Ca2+ fluxes, and sarcoplasmic reticulum loading in pressure-overload hypertrophy
    1. L Delbridge
    2. H Satoh
    3. W Yuan
    4. J Bassani
    5. M Qi
    6. KS Ginsburg
    7. AM Samarel
    8. DM Bers
    (1997)
    American Journal of Physiology-Heart and Circulatory Physiology 272:H2425–H2435.
  12. 12
  13. 13
    Ionic mechanism of electrical alternans
    1. JJ Fox
    2. JL McHarg
    3. RF Gilmour
    (2002)
    American Journal of Physiology-Heart and Circulatory Physiology 282:H516–H530.
    https://doi.org/10.1152/ajpheart.00612.2001
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63

Decision letter

  1. Richard Aldrich
    Senior Editor; The University of Texas at Austin, United States
  2. José D Faraldo-Gómez
    Reviewing Editor; National Heart, Lung and Blood Institute, National Institutes of Health, United States
  3. Blanca Rodriguez
    Reviewer; Department of Computer Science, University of Oxford, Oxford, United Kingdom
  4. Colleen Clancy
    Reviewer

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Variability and compensation of cardiomyocyte ionic conductances at the population level" for consideration by eLife. Your article has been reviewed by Richard Aldrich as the Senior Editor, a Reviewing Editor, and three reviewers. The following individuals involved in review of your submission have agreed to reveal their identity: David Christini (Reviewer #1), Blanca Rodriguez (Reviewer #2); Colleen Clancy (Reviewer #3).

The reviewers' comments are appended below. As you can see, the reviewers raise a number of questions and concerns pertaining the methodology used in your study and the interpretation and presentation of your results. Based on this assessment, we cannot accept your manuscript for publication in its present form. Nonetheless, the reviewers and the Reviewing Editor concur that the issues raised can be resolved in a reasonable timeframe, and therefore we would like to encourage to resubmit a revised version of your manuscript that addresses the reviewers' comments, along with a point-by-point response

Reviewer #1:

The authors present an intriguing study that adds new insight to a developing area of research utilizing rich electrophysiological and/or clinical data sets to reparameterize cellular cardiac models.

In this study, good enough solution (GES) parameterizations are first developed in simulations in which calcium dynamic properties were constrained by requiring minimization of an error function of calcium and sodium concentration differences. Such constraints produced model variants that illuminated a very interesting compensatory relationship between inward Ca2+ and outward K+ currents. These findings were used as evidence in support of the "central hypothesis that parameters are constrained predominantly by features of the CaT". I am confused by the logic flow that led to that hypothesis. One issue that comes to mind is that in real cells, how can the authors know if it is CaT that is controlling (via feedback) the ionic current expression (the suggested mechanism), OR if it is the other way around, e.g., inward Ca2+ channel expression governing outward K+ channel expression (and/or vice versa, via any number of unknown second-messenger mechanisms), which might result in consistent (constrained) CaT?

In several locations, especially in the background and conclusions, the authors emphasize that by focusing on CaT as the main model constraining variables, this study presents a fundamentally new approach. In my view, the approach used in this study is a new twist, rather than a fundamental advancement. Although I did not reread the papers (cited in the background) that used related approaches, my recollection is that some of them at least refer to the idea that adding calcium measurements to their approaches would add value. I thus see this approach as a straightforward variant of previous approaches.

By focusing so much attention on the newness of the CaT constraint, I think the authors distract from two aspects that are more likely to be of interest to the audience: (1) the compensatory Ca2+ / K+ findings and (2) the utilization of the multiple species of the Hybrid Mouse Diversity Panel. Both of these are strengths of the paper.

It is not clear to me why the authors choose not to include membrane potential in the error function. As they note in subsection “Computationally determined good enough solutions”, "the AP waveforms exhibit larger variations owing to the fact that the GES search does not involve any voltage sensing." Because of this choice, AP waveforms are often non-physiological (e.g. the huge variation shown in Figure 3B), which is a notable weakness.

Overall, this study is intriguing, notwithstanding the aforementioned issues.

Reviewer #2:

The primary aim of the paper is to evaluate how parameters in a mouse cardiac electrophysiology model would be constrained using calcium transient values rather than action potential values. They show that, in mice, constraining the model search using calcium transient yields a correlation between the sum of two potassium currents and the main calcium current.

The idea is interesting. The paper however is superficial in the results provided and also the discussion. This is primarily related to the fact that a comparison of results with action potential calibration versus calcium calibration is not shown, and the results seem rather species dependent (to mice). Both aspects could be easily explored in more depth. I also have further suggestions to better position this work in the context of previous relevant studies in cardiac electrophysiology.

The paper postulates that "model parameters are predominantly constrained by feedback sensing of Ca2+, and potentially other ions (e.g. Na+) affecting ion channel regulation".

The most interesting result in the paper is the correlation between the sum of two potassium currents and ICa,L, demonstrated both computationally and experimentally. Even though correlations between 2 and 3 parameters have been previously explored and this could be mentioned in the paper (e.g. Sanchez et al., 2014, Britton et al., 2013, Muszkiewicz et al., 2018), their sum hasn't been considered and this seems important. I have two general comments about this:

1) Only 6 parameters are varied and, of these, only 3 conductances are shown to be correlated. The statement above therefore only applies to a subset of conductances. Moreover, previous work on populations of models has shown that certain conductances are indeed constrained when calibrating only with voltage biomarkers (for example the sodium current by the upstroke velocity). How does this fit in the primary hypothesis of the paper? Would those conductances have voltage sensors? Please add discussion on this.

2) The analysis of the AP phenotypes following the Ca-based calibration would be important (it is in fact mentioned in the Discussion but not shown in the results). A previous study shows that AP calibration leads to variable calcium phenotypes in human atrial cells (Muszkiewicz et al., 2018). Something similar can happen with Ca-based calibration. The analysis would be both feasible and very interesting.

3) What does AP calibration (in addition to Calcium calibration) yield in terms of ionic conductances correlation and constrain? The study does not demonstrate that AP calibration doesn't achieve the same sort of correlation in the ionic conductances than calcium transient calibration. This would be important to demonstrate the main premise of the paper.

4) This paper focuses on mice (rather than human or other species). This should certainly be highlighted in the title and somewhat detracts from the relevance of the study. The repolarization ionic currents, action potential and calcium waveforms are very different in mice and human, particular from the ventricles. This is an aspect that should be analysed and discussed in the paper.

It would indeed be interesting to test whether a similar level of conductances correlation is observed or not when using human models (for atria and ventricles) calibrated with calcium data. How many repolarization currents would need to be considered in the sum to obtain correlation with the calcium current? This is feasible given the availability of human models and calcium data, and the fact that the models have already been used for population studies on different human cell types (ventricular, atrial and iPS, see references mentioned above). A similar comparison was for example done in the study by Devenyi and Sobie, 2016.

5) I am confused by what is shown in Figure 5, and how it contributes to the aim of the paper. It shows that electrotonic coupling reduces APD heterogeneity (as shown previously in the context of alternans Zaniboni et al., 2016 as well as with populations of models for example in Sanchez et al., 2017). The paper states that calcium amplitude heterogeneity is also reduced (subsection “Compensation at the organ scale”). However, I see a very small reduction of calcium amplitude heterogeneity in Figure 5E and I certainly don't see a very strong effect. I would expect however spatial heterogeneity of calcium amplitude to be reduced by calcium flow through gap junctions. This model feature seems crucial in this context. In addition, APD calibration would also achieve the same outputs, wouldn't it?

6) The section on "Cardiac hypertrophic response to a stressor" contains the interesting correlation shown in Figure 6 and then it is mostly speculative. I don't see a clear contribution towards the goals of the paper as set in the introduction.

Reviewer #3:

In paper from Rees and colleagues, the authors present a model to predict variability of ion channel conductances in a population of subjects with variable cardiac gene expression patterns. Unlike previous models that account for variability using features of the action potential, the new model focuses on the calcium transient (CaT).

The authors build on previous work demonstrating that many distinct parameter sets can elicit the same phenomenological behavior. Specifically, the authors point out that most previous works have focused on attributes of transmembrane potential as the important phenomenological behaviors to be replicated by all parameter sets. The authors propose that sensory mechanisms exist for monitoring the calcium transient and intracellular sodium dynamics. The authors vary six ionic conductances and examine what combinations replicate the calcium transient and intracellular sodium dynamics of the baseline model. Through this analysis, they are able to identify parameters that are highly constrained, relatively insensitive, or highly correlated with other parameters. As pointed out by the authors, this novel work using the calcium transient and intracellular sodium dynamics as opposed to transmembrane potential dynamics to constrain the model is important, because when models are fit to replicate physiological phenomena, the final model parameters critically depend on what physiological characteristics are used to constrain the model. To test their model, the authors used experiments using the Hybrid Mouse Diversity Panel.

The authors state at the end of their discussion, that this work has very broad implications in the realm of understanding how drugs will affect cardiac health in a diverse population. Since ion channel expression levels can span a wide range (both among individuals in a population as well as between cells in a single individual), the effect of a drug that binds to a specific ion channel could vary widely between individuals or from cell to cell. Therefore, to accurately characterize the effects of a drug that interacts with cardiac ion channels across an entire population, it is important to fully characterize the parameter space of physiologically allowable ionic conductances. This is an important point and could be emphasized more.

What seems to be missing from the study at present is a description of the fundamental insight that is gained by the Ca-based GES search compared to the V-based GES search. It's hard to imagine that the V-based search wouldn't reveal the relationship between Ito,f and IKur and ICa,L. Although I understand that the Ca signal is (likely) more directly related to ion-channel regulation, there is no particular control mechanism included in the analysis – so what does the current state analysis provide that the V-based doesn't? It would be so useful to show a direct comparison of the Ca- and V-based GES searches. It should also be mentioned that there are a host of regulatory mechanisms not accounted for in either search method.

1) More discussion of the range of APs waveforms observed would be helpful, since that is what has been examined in all other cardiac AP population approaches. Additionally, in the Discussion section ("physiological AP waveform whose duration matches the experimentally observed variation in isolated myocytes even though the voltage signal is not used to constrain model parameters"), doesn't seem to be entirely supported by the data shown. Figure 3B would better support this argument if the full range of simulated AP behaviors were quantified and compared to experimental AP data. It is not clear how the APs shown in Figure 3B (and later in Figure 5) were chosen, or if they represent the full range of behaviors observed with this approach.

2) Similarly, for Figure 5C, it would be clearer if cell-to-cell variability in the simulated APs was quantified and compared that to the experimental range. It seems that the variability in the AP waveforms shown in Figure 5C represent a much wider range of morphologies (in APD90, for example), than in the experimental APs shown in Figure 5—figure supplement 1.

3) It is not clear what is defined as the "normal CaT (Calcium Transient)", or how that was decided. The reference values (Gref) are defined as the conductance in the baseline model, which produce the normal CaT. Is the normal CaT the same as the reference outputs defined in Table 6? Are these the only two outputs used to define the normal CaT? Also, for the reference CaT outputs (as shown in Table 6), how were these values decided? Include a citation (or justification) for the source of the reference/normal values.

4) In the populations shown in Figure 5, is a similar coupling effect seen for average cytosolic Ca2+ over one beat? Does this range of simulated values, for either CaT measure, mimic the experimental range?

5) My understanding is that the point of Figure 5 is that even if the conductances considered here aren't tightly correlated, the Ca2+ transient and especially the AP are constrained due to electrical coupling in tissue. However, what is not clear to me based on the authors' discussion is if electrical coupling is enough to constrain the Ca2+ transient and AP to within a healthy range. Based on looking at the figure, electrical coupling does not constrain the Ca2+ transient to within the 5% error used in the earlier good enough solutions. But how well constrained do the Ca2+ transient and AP really need to be? Also, how would things change if the parameter sets of the cells in the tissue came from the good enough solution?

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Variability and compensation of cardiomyocyte ionic conductances at the population level" for further consideration at eLife. Your revised article has been favorably evaluated by Richard Aldrich (Senior Editor), José D. Faraldo-Gómez as Reviewing Editor, and three reviewers.

The manuscript has been improved but there are some remaining issues that we would like you to address before acceptance, as outlined below. Although your revised manuscript will not be sent to the reviewers, we encourage you to carefully consider their suggestions.

Reviewer #1:

The authors have done a good job of addressing the general concerns that I raised in my original review. Notably, they have sharpened up the presentation of the overall hypothesis of CaT regulation, via blood-pressure feedback, serving as the mechanism governing underlying levels of Ca-cycling proteins and the action potential. As I noted in my original review, the computational findings are supported by analyses of the Hybrid Mouse Diversity Panel, which strengthen the study significantly and set it apart from a purely computational study.

One concern I have relates to a point made in response to my review, as well as the other reviewers (and in the Discussion section of the revised manuscript). The authors provide results and analyses to defend the approach of focusing on calcium for the GES determination, rather than also incorporating the action potential. They argue that the variation in the action potential waveform that occurs using the CaT-based parameterization is within expected physiological variability, especially if we expand our definition of physiological variability to include outlier AP waveforms that are excluded via selection bias by experimentalists who "often" discard sick-looking AP waveforms. My first problem here is that this statement is an over simplification of what good experimentalists do and why. Sure, some may just exclude APs that look wrong, but that is not good accepted practice. Instead, they use specific quantifiable criteria, such as minimum resting membrane potential, which have been demonstrated to be correlated to damaged cells, for exclusion. While such criteria may still be excluding physiological APs, my point is that it is more nuanced than the text implies and, I think, does a disservice to the actual experiments and the reported data. Beyond that issue, within the very same paragraph, the authors note that the variability that occurs in CaT using voltage-based parameterization is "highly variable and often nonphysiological." So why is it that the variability that arises using the approach that fits with the authors hypothesis is physiological, whereas the variability using the approach that does not fit the authors hypothesis is "nonphysiological?" Put another way, when the authors define nonphysiological CaT, are they assuming that experimentalists who recorded and reported CaT have never thrown away "sick looking" calcium transient data, which, if included in the literature, would expand the range of physiological? This point is really just about framing, rather than about the overall justification of the approach, as I do agree with the overall conclusion that these findings support using CaT. But I think that point can be made without the "sick looking" discussion or expanding that discussion with more context.

I agree with the point made by referee #2 that it would be useful for the authors to deposit the code used to generate the results.

Reviewer #2:

The authors have followed most of our suggestions and the paper is stronger. The title remains too generic and would benefit from at least mentioning calcium as the main novelty of the paper (and perhaps mice too). Otherwise it could very well be another review on the topic.

https://doi.org/10.7554/eLife.36717.025

Author response

Reviewer #1:

The authors present an intriguing study that adds new insight to a developing area of research utilizing rich electrophysiological and/or clinical data sets to reparameterize cellular cardiac models.

In this study, good enough solution (GES) parameterizations are first developed in simulations in which calcium dynamic properties were constrained by requiring minimization of an error function of calcium and sodium concentration differences. Such constraints produced model variants that illuminated a very interesting compensatory relationship between inward Ca2+ and outward K+ currents. These findings were used as evidence in support of the "central hypothesis that parameters are constrained predominantly by features of the CaT". I am confused by the logic flow that led to that hypothesis. One issue that comes to mind is that in real cells, how can the authors know if it is CaT that is controlling (via feedback) the ionic current expression (the suggested mechanism), OR if it is the other way around, e.g., inward Ca2+ channel expression governing outward K+ channel expression (and/or vice versa, via any number of unknown second-messenger mechanisms), which might result in consistent (constrained) CaT?

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 Ca current-K current ratio. We have revised the text in the introduction to more clearly state this hypothesis.

In several locations, especially in the background and conclusions, the authors emphasize that by focusing on CaT as the main model constraining variables, this study presents a fundamentally new approach. In my view, the approach used in this study is a new twist, rather than a fundamental advancement. Although I did not reread the papers (cited in the background) that used related approaches, my recollection is that some of them at least refer to the idea that adding calcium measurements to their approaches would add value. I thus see this approach as a straightforward variant of previous approaches.

We believe that the physiological hypothesis, as articulated above, is the novel insight from this study, rather than the technical approach underlying the GES search, which we acknowledge uses standard methods. We have made this clear in the text.

By focusing so much attention on the newness of the CaT constraint, I think the authors distract from two aspects that are more likely to be of interest to the audience: (1) the compensatory Ca2+ / K+ findings and (2) the utilization of the multiple species of the Hybrid Mouse Diversity Panel. Both of these are strengths of the paper.

It is not clear to me why the authors choose not to include membrane potential in the error function. As they note in subsection “Computationally determined good enough solutions”, "the AP waveforms exhibit larger variations owing to the fact that the GES search does not involve any voltage sensing." Because of this choice, AP waveforms are often non-physiological (e.g. the huge variation shown in Figure 3B), which is a notable weakness.

Overall, this study is intriguing, notwithstanding the aforementioned issues.

We have given more emphasis to the compensatory Ca/K findings and the utilization of the Hybrid Mouse Diversity Panel in the introduction. We note that although AP waveforms are uniform in tissue, there is a marked variation in isolated ventricular myocytes. Moreover, experimentalists often exclude atypical “sick-looking” AP waveforms in patch clamp experiments, which are assumed to be “non-physiological” due to cell damage during enzymatic isolation or patching. Thus, AP waveform variation in single myocytes may be underestimated because of this selection bias. In coupled tissue, however, occasional atypical AP waveforms would be voltage-clamped by their “healthy” AP waveform neighbors to produce an overall uniform AP waveform. If we had included voltage sensors together with Ca2+ sensors, the AP waveforms would, like the CaT, have been much more uniform. However, we did not feel it necessary, since (1) the physiological basis for voltage sensors is not clear, whereas the physiological basis for the CaT as a surrogate for blood pressure mediated feedback to the heart is well-known; (2) if we had included voltage sensors to constrain AP waveform, its variation could become arbitrarily narrower than the observed AP waveform variation in isolated myocytes. In response to this comment and to the similar first major comment of Reviewer #2, we have discussed these points and added new simulation results in which the Ca2+ sensors were substituted for voltage sensors (new Figure 3—figure supplement 3). Although the AP waveform was readily constrained, the CaT was highly variable and often nonphysiological, and the Ca/K current ratio relationship was no longer preserved, since other inward and outward currents (in particular the Na+-Ca2+ exchanger) could regulate AP duration when the CaT was not constrained.

Reviewer #2:

The primary aim of the paper is to evaluate how parameters in a mouse cardiac electrophysiology model would be constrained using calcium transient values rather than action potential values. They show that, in mice, constraining the model search using calcium transient yields a correlation between the sum of two potassium currents and the main calcium current.

The idea is interesting. The paper however is superficial in the results provided and also the discussion. This is primarily related to the fact that a comparison of results with action potential calibration versus calcium calibration is not shown, and the results seem rather species dependent (to mice). Both aspects could be easily explored in more depth. I also have further suggestions to better position this work in the context of previous relevant studies in cardiac electrophysiology.

As already mentioned in the last response to reviewer #1, we have added new simulation results voltage sensing alone (new Figure 3—figure supplement 3) that can be readily compared to Ca2+ sensing alone. With voltage sensing alone, the AP waveform was readily constrained as expected, but the CaT was highly variable and often nonphysiological. The Ca/K current ratio relationship was no longer preserved, since other inward and outward currents could regulate AP duration when the CaT was not constrained. As noted above, our hypothesis provides a clear physiological link for Ca2+ sensing, 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. This provides a very 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 Ca current-K current ratio. 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 now shown in the new results. We have emphasized these points in the revised Discussion section.

The paper postulates that "model parameters are predominantly constrained by feedback sensing of Ca2+, and potentially other ions (e.g. Na+) affecting ion channel regulation". The most interesting result in the paper is the correlation between the sum of two potassium currents and ICa,L, demonstrated both computationally and experimentally. Even though correlations between 2 and 3 parameters have been previously explored and this could be mentioned in the paper (e.g. Sanchez et al., 2014, Britton et al., 20134, Muszkiewicz et al., 2018), their sum hasn't been considered and this seems important. I have two general comments about this:

We have mentioned the additional references in the revised paper.

1) Only 6 parameters are varied and, of these, only 3 conductances are shown to be correlated. The statement above therefore only applies to a subset of conductances. Moreover, previous work on populations of models has shown that certain conductances are indeed constrained when calibrating only with voltage biomarkers (for example the sodium current by the upstroke velocity). How does this fit in the primary hypothesis of the paper? Would those conductances have voltage sensors? Please add discussion on this.

The 6 parameters that we used play a major role in regulating the CaT, based on our hypothesis that the CaT is a surrogate for blood pressure that feeds back to the heart via the autonomic nervous system. The reviewer makes a good point about the Na current, since 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. Thus, it seems reasonable to propose that Na current density (and gap junction coupling) could also be regulated by Ca2+ sensing to ensure that conduction velocity is adequate to generate a synchronous blood pressure waveform. We have added these points to the Discussion section.

2) The analysis of the AP phenotypes following the Ca-based calibration would be important (it is in fact mentioned in the Discussion but not shown in the results). A previous study shows that AP calibration leads to variable calcium phenotypes in human atrial cells (Muszkiewicz et al., 2018). Something similar can happen with Ca-based calibration. The analysis would be both feasible and very interesting.

Our results show that CaT calibration results in considerable AP waveform variability (Figure 5C) as expected if the AP waveform is not constrained. While this is consistent with the converse finding by Muszkiewicz et al., that the CaT is highly variable when the AP waveform is constrained, constraining the AP waveform does not produce the correlation between Ca and K currents observed in the HMDP, and is difficult to justify biologically given the absence of voltage sensing mechanism. At the tissue level, AP waveform variability is greatly smoothed by diffusive coupling, and so does not require tight control at the single myocyte level. As discussed above, the single myocyte AP variability is considerable when measured experimentally in patch clamp studies. We have added those points to the Discussion section.

3) What does AP calibration (in addition to Calcium calibration) yield in terms of ionic conductances correlation and constrain? The study does not demonstrate that AP calibration doesn't achieve the same sort of correlation in the ionic conductances than calcium transient calibration. This would be important to demonstrate the main premise of the paper.

As noted above, we have added new simulations directly comparing Ca2+ sensing alone with voltage sensing alone (Figure 3—figure supplement 3). With voltage sensing alone, the AP waveform was readily constrained as expected, but the CaT was highly variable and often nonphysiological. The Ca/K current ratio relationship was no longer preserved, since other inward and outward currents could regulate AP duration when the CaT was not constrained.

4) This paper focuses on mice (rather than human or other species). This should certainly be highlighted in the title and somewhat detracts from the relevance of the study. The repolarization ionic currents, action potential and calcium waveforms are very different in mice and human, particular from the ventricles. This is an aspect that should be analysed and discussed in the paper.

It would indeed be interesting to test whether a similar level of conductances correlation is observed or not when using human models (for atria and ventricles) calibrated with calcium data. How many repolarization currents would need to be considered in the sum to obtain correlation with the calcium current? This is feasible given the availability of human models and calcium data, and the fact that the models have already been used for population studies on different human cell types (ventricular, atrial and iPS, see references mentioned above). A similar comparison was for example done in the study by Devenyi and Sobie, 2016.

We have performed a similar GES study using a rabbit ventricular myocyte model, in which the AP and EC coupling features are closer to human. Although we obtained similar results, we elected not to include these findings in the current manuscript to avoid making it unduly long, and also because a rabbit analog of the HMDP is not available to validate the computational results at a level comparable to what has been achieved for mice here. However, we plan to submit the rabbit model results in a separate manuscript. The results of this GES search is illustrated above for the reviewers’ benefit. We have also added a paragraph in the discussion to suggest the use of cardiomyocytes (CMs) derived from induced pluripotent stem cells (iPSC) as an alternative to the HMDP to extend the present study to human. 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, if present and statistically distinguishable between 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.

5) I am confused by what is shown in Figure 5, and how it contributes to the aim of the paper. It shows that electrotonic coupling reduces APD heterogeneity (as shown previously in the context of alternans Zaniboni et al., 2016 as well as with populations of models for example in Sanchez et al., 2017). The paper states that calcium amplitude heterogeneity is also reduced (subsection “Compensation at the organ scale”). However, I see a very small reduction of calcium amplitude heterogeneity in Figure 5E and I certainly don't see a very strong effect. I would expect however spatial heterogeneity of calcium amplitude to be reduced by calcium flow through gap junctions. This model feature seems crucial in this context. In addition, APD calibration would also achieve the same outputs, wouldn't it?

We agree that unlike APD heterogeneity, electrotonic coupling had a more modest effect on CaT heterogeneity, although it was still significant. It is possible that the heterogeneity would be further reduced by Ca flow through gap junctions, but we did not simulate this possibility. We doubt that this effect would be very significant since Ca diffusion through gap junctions is orders of magnitude slower than voltage diffusion.

6) The section on "Cardiac hypertrophic response to a stressor" contains the interesting correlation shown in Figure 6 and then it is mostly speculative. I don't see a clear contribution towards the goals of the paper as set in the introduction.

We acknowledge that this section is very speculative. However, we prefer to keep it as an example of how the approach outlined in the paper could be applied to obtain insights (or stimulate future avenues of investigation) relevant to disease conditions.

Reviewer #3:

In paper from Rees and colleagues, the authors present a model to predict variability of ion channel conductances in a population of subjects with variable cardiac gene expression patterns. Unlike previous models that account for variability using features of the action potential, the new model focuses on the calcium transient (CaT).

The authors build on previous work demonstrating that many distinct parameter sets can elicit the same phenomenological behavior. Specifically, the authors point out that most previous works have focused on attributes of transmembrane potential as the important phenomenological behaviors to be replicated by all parameter sets. The authors propose that sensory mechanisms exist for monitoring the calcium transient and intracellular sodium dynamics. The authors vary six ionic conductances and examine what combinations replicate the calcium transient and intracellular sodium dynamics of the baseline model. Through this analysis, they are able to identify parameters that are highly constrained, relatively insensitive, or highly correlated with other parameters. As pointed out by the authors, this novel work using the calcium transient and intracellular sodium dynamics as opposed to transmembrane potential dynamics to constrain the model is important, because when models are fit to replicate physiological phenomena, the final model parameters critically depend on what physiological characteristics are used to constrain the model. To test their model, the authors used experiments using the Hybrid Mouse Diversity Panel.

The authors state at the end of their discussion, that this work has very broad implications in the realm of understanding how drugs will affect cardiac health in a diverse population. Since ion channel expression levels can span a wide range (both among individuals in a population as well as between cells in a single individual), the effect of a drug that binds to a specific ion channel could vary widely between individuals or from cell to cell. Therefore, to accurately characterize the effects of a drug that interacts with cardiac ion channels across an entire population, it is important to fully characterize the parameter space of physiologically allowable ionic conductances. This is an important point and could be emphasized more.

What seems to be missing from the study at present is a description of the fundamental insight that is gained by the Ca-based GES search compared to the V-based GES search. It's hard to imagine that the V-based search wouldn't reveal the relationship between I_to,f+I_Kur and I_CaL. Although I understand that the Ca signal is (likely) more directly related to ion-channel regulation, there is no particular control mechanism included in the analysis – so what does the current state analysis provide that the V-based doesn't? It would be so useful to show a direct comparison of the Ca- and V-based GES searches. It should also be mentioned that there are a host of regulatory mechanisms not accounted for in either search method.

As noted in the response to reviewer 2, we have added new simulations directly comparing Ca2+ sensing alone with voltage sensing alone (Figure 3—figure supplement 3). With voltage sensing alone, the AP waveform was readily constrained as expected, but the CaT was highly variable and often nonphysiological. The Ca/K current ratio relationship was no longer preserved, since other inward and outward currents could regulate AP duration when the CaT was not constrained. As noted above, our hypothesis provides a clear physiological link for Ca2+ sensing, 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. This provides a very 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 Ca current-K current ratio. 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 now shown in the new results. We have emphasized these points in the revised text.

1) More discussion of the range of APs waveforms observed would be helpful, since that is what has been examined in all other cardiac AP population approaches. Additionally, in the Discussion section ("physiological AP waveform whose duration matches the experimentally observed variation in isolated myocytes even though the voltage signal is not used to constrain model parameters"), doesn't seem to be entirely supported by the data shown. Figure 3B would better support this argument if the full range of simulated AP behaviors were quantified and compared to experimental AP data. It is not clear how the APs shown in Figure 3B (and later in Figure 5) were chosen, or if they represent the full range of behaviors observed with this approach.

As noted above in the response from reviewer 2, it is clear from visual inspection of the experimental traces in Figure 5—figure supplement 1 that there is marked variation in the AP waveform among different myocytes, ranging from about 10 to 80 ms. This is comparable to the range of APD values in Figure 5E, although there is not sufficient experimental data to construct a histogram for direct comparison. In the text, we have modified the text to state: “A remarkable and nontrivial finding of the present computational study is that Ca2+ 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.”

2) Similarly, for Figure 5C, it would be clearer if cell-to-cell variability in the simulated APs was quantified and compared that to the experimental range. It seems that the variability in the AP waveforms shown in Figure 5C represent a much wider range of morphologies (in APD90, for example), than in the experimental APs shown in Figure 5—figure supplement 1.

See the response to the previous comment.

3) It is not clear what is defined as the "normal CaT (Calcium Transient)", or how that was decided. The reference values (Gref) are defined as the conductance in the baseline model, which produce the normal CaT. Is the normal CaT the same as the reference outputs defined in Table 6? Are these the only two outputs used to define the normal CaT? Also, for the reference CaT outputs (as shown in Table 6), how were these values decided? Include a citation (or justification) for the source of the reference/normal values.

We arbitrarily constrained the steady-state CaT amplitude, [Ca]i, time averaged

cytosolic Ca2+concentration, and the intracellular sodium concentration [Na]ito fall within a 5% tolerance of the control values at a 4 Hz pacing frequency. Slightly increasing or decreasing the tolerance only increases or decreases the range of conductances producing a normal CaT but does not destroy the correlation between Ca and K currents shown in Figures 3F and Figure 4C. Furthermore, as explained in response to a similar question of reviewer #2, we chose a 4 Hz pacing frequency based on the fact that the Bondarenko et al., ventricular myocyte mouse model used in the present study was benchmarked against experimental measurements of CaT amplitude (Figure 18 in Bondarenko et al., 2018) for frequencies ranging from 0.5 to 6 Hz. Since different model parameters corresponding to different strains reproduce similar CaT amplitudefrequency curves as the Bondarenko et al., model, we do not expect the choice of pacing frequency to be critically important for calibrating model parameters that produce a normal electrophysiological phenotype. We have discussed this point in the Discussion section.

4) In the populations shown in Figure 5, is a similar coupling effect seen for average cytosolic Ca2+ over one beat? Does this range of simulated values, for either CaT measure, mimic the experimental range?

Yes, a similar coupling effect is seen for the average cytosolic Ca2+ concentration and is stronger than for the CaT amplitude. We have added a new Figure (Figure 5—figure supplement 2) to show this effect. We did not measure CaTs experimentally, which are very hard to quantify accurately to compare with the simulated data.

5) My understanding is that the point of Figure 5 is that even if the conductances considered here aren't tightly correlated, the Ca2+ transient and especially the AP are constrained due to electrical coupling in tissue. However, what is not clear to me based on the authors' discussion is if electrical coupling is enough to constrain the Ca2+ transient and AP to within a healthy range. Based on looking at the figure, electrical coupling does not constrain the Ca2+ transient to within the 5% error used in the earlier good enough solutions. But how well constrained do the Ca2+ transient and AP really need to be? Also, how would things change if the parameter sets of the cells in the tissue came from the good enough solution?

We apologize for the confusion here. In the GES search, the CaT was constrained to within 5%. However, in the patch clamp experiments, the K and Ca currents had to be measured in different myocytes to isolate the currents properly, so that we could not verify whether the Ca/K current ratio was constrained at the single myocyte level, even though the average ratio for all of the myocytes studied was constant. In this simulation, we addressed the issue of whether a constant Ca/K current ratio at the single myocyte level was important or not when extrapolated to the tissue level. Therefore, we assigned random experimentally measured values of the Ca current and K current densities for the C57BL/6 strain (blue point cluster in Figure 5B) to the AP model to generate the AP and CaT distributions shown in Figure 5E. Because the Ca/K current ratio was not constrained in this population of models, both APD and CaT showed marked variation (unlike the case for the GES search in which the CaT was within 5%). Under these conditions, the APD variance was markedly reduced by coupling, and the CaT variance was modestly reduced. Thus, this tissue simulation represents a worse-case scenario in which the Ca/K current ratio was unconstrained at the single myocyte level, but the mean Ca/K current ratio from many myocytes was constrained. The conclusion is that in coupled tissue, marked variation is the APD can be compensated for at the tissue level by coupling, and marked variation in the CaT can be somewhat compensated for by coupling. However, although we could not ascertain whether the Ca/K current ratio is constant at the single myocyte level to produce low CaT variance (e.g. within 5%), the cell shortening data in Figure 4B shows that in the C57BL/6 strain, cell shortening (as a surrogate for CaT amplitude) varied about 2-fold. This is considerably less than the variance in the CaT amplitude distribution for the uncoupled simulated AP models in Figure 5E, in which the Ca/K current ratio was completely unconstrained. We have revised subsection “Compensation at the organ scale” to clarify these points.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that we would like you to address before acceptance, as outlined below. Although your revised manuscript will not be sent to the reviewers, we encourage you to carefully consider their suggestions.

Reviewer #1:

The authors have done a good job of addressing the general concerns that I raised in my original review. Notably, they have sharpened up the presentation of the overall hypothesis of CaT regulation, via blood-pressure feedback, serving as the mechanism governing underlying levels of Ca-cycling proteins and the action potential. As I noted in my original review, the computational findings are supported by analyses of the Hybrid Mouse Diversity Panel, which strengthen the study significantly and set it apart from a purely computational study.

One concern I have relates to a point made in response to my review, as well as the other reviewers (and in the Discussion section of the revised manuscript). The authors provide results and analyses to defend the approach of focusing on calcium for the GES determination, rather than also incorporating the action potential. They argue that the variation in the action potential waveform that occurs using the CaT-based parameterization is within expected physiological variability, especially if we expand our definition of physiological variability to include outlier AP waveforms that are excluded via selection bias by experimentalists who "often" discard sick-looking AP waveforms. My first problem here is that this statement is an over simplification of what good experimentalists do and why. Sure, some may just exclude APs that look wrong, but that is not good accepted practice. Instead, they use specific quantifiable criteria, such as minimum resting membrane potential, which have been demonstrated to be correlated to damaged cells, for exclusion. While such criteria may still be excluding physiological APs, my point is that it is more nuanced than the text implies and, I think, does a disservice to the actual experiments and the reported data. Beyond that issue, within the very same paragraph, the authors note that the variability that occurs in CaT using voltage-based parameterization is "highly variable and often nonphysiological." So why is it that the variability that arises using the approach that fits with the authors hypothesis is physiological, whereas the variability using the approach that does not fit the authors hypothesis is "nonphysiological?" Put another way, when the authors define nonphysiological CaT, are they assuming that experimentalists who recorded and reported CaT have never thrown away "sick looking" calcium transient data, which, if included in the literature, would expand the range of physiological? This point is really just about framing, rather than about the overall justification of the approach, as I do agree with the overall conclusion that these findings support using CaT. But I think that point can be made without the "sick looking" discussion or expanding that discussion with more context.

We agree with the reviewer that our main point that the CaT suffices to regulate ion channel conductances without voltage sensing can be made without invoking “sick looking” cells. Accordingly, we have replaced the part of the text in the Discussion section:

“Experimentalists often exclude atypical “sick-looking'' AP waveforms in patch clamp experiments, which are assumed to be “non-physiological'' due to cell damage during enzymatic isolation or patching. Thus, AP waveform variation in isolated myocytes may be underestimated because of this selection bias. In coupled tissue, however, occasional atypical AP waveforms would be voltage-clamped by their “healthy'' AP waveform neighbors to produce an overall uniform AP waveform.”

by

“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 is greatly reduced in tissue because less frequent atypical AP waveforms are voltage-clamped by the more typical AP waveforms of their neighbors.”

I agree with the point made by referee #2 that it would be useful for the authors to deposit the code used to generate the results.

The codes that produced the results have been uploaded to a public repository cited in the article as

Code availability

The codes used in this work are available at:

https://github.com/circs/GES

Reviewer #2:

The authors have followed most of our suggestions and the paper is stronger. The title remains too generic and would benefit from at least mentioning calcium as the main novelty of the paper (and perhaps mice too). Otherwise it could very well be another review on the topic.

Following the suggestion of the reviewer, we have modified the title to refer explicitly to the calcium transient as the main novelty of our study and to specify that it was carried out in mouse populations (i.e. using the Hybrid Mouse Diversity Panel).

https://doi.org/10.7554/eLife.36717.026

Article and author information

Author details

  1. Colin M Rees

    1. Physics Department, Northeastern University, Boston, United states
    2. Center for Interdisciplinary Research on Complex Systems, Northeastern University, Boston, United States
    Contribution
    Software, Formal analysis, Investigation
    Contributed equally with
    Jun-Hai Yang
    Competing interests
    No competing interests declared
  2. Jun-Hai Yang

    1. Department of Medicine (Cardiology), Cardiovascular Research Laboratory, David Geffen School of Medicine, University of California, Los Angeles, United states
    2. Department of Physiology, David Geffen School of Medicine, University of California, Los Angeles, United States
    Contribution
    Data curation, Investigation, Writing—review and editing, Made a very substantial contribution to this work, Performed all the patch clamp experiments and contractility measurements, Recorded three major ionic currents and cell shortening for a large number of cells and different mice strains, Provided key data to validate computational modeling predictions, Responsible for curation of the data, Helped review the draft of the paper
    Contributed equally with
    Colin M Rees
    Competing interests
    No competing interests declared
  3. Marc Santolini

    1. Physics Department, Northeastern University, Boston, United states
    2. Center for Interdisciplinary Research on Complex Systems, Northeastern University, Boston, United States
    Contribution
    Data curation, Validation, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1491-0120
  4. Aldons J Lusis

    1. Department of Medicine (Cardiology), Cardiovascular Research Laboratory, David Geffen School of Medicine, University of California, Los Angeles, United states
    2. Department of Physiology, David Geffen School of Medicine, University of California, Los Angeles, United States
    3. Department of Microbiology, David Geffen School of Medicine, University of California, Los Angeles, United States
    4. Department of Human Genetics, David Geffen School of Medicine, University of California, Los Angeles, United States
    Contribution
    Resources, Funding acquisition
    Competing interests
    No competing interests declared
  5. James N Weiss

    1. Department of Medicine (Cardiology), Cardiovascular Research Laboratory, David Geffen School of Medicine, University of California, Los Angeles, United states
    2. Department of Physiology, David Geffen School of Medicine, University of California, Los Angeles, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition
    Competing interests
    No competing interests declared
  6. Alain Karma

    1. Physics Department, Northeastern University, Boston, United states
    2. Center for Interdisciplinary Research on Complex Systems, Northeastern University, Boston, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Methodology, Writing—original draft
    For correspondence
    a.karma@northeastern.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7032-9862

Funding

National Heart, Lung, and Blood Institute (5R01HL114437)

  • Colin M Rees
  • Jun-Hai Yang
  • Marc Santolini
  • Aldons J Lusis
  • James N Weiss
  • Alain Karma

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This research was supported by NIH/NHLBI grant 5R01HL114437 and by the Laubisch and Kawata Endowments. The authors thank Yibin Wang for valuable input into this study. AK acknowledges stimulating discussions with Eve Marder in the early stage of this work.

Ethics

Animal experimentation: This study was approved by the UCLA Chancellor's Animal Research Committee (ARC 2003-063-23B) and performed in accordance with the Guide for the Care and Use of Laboratory Animals published by the United States National Institutes of Health (NIH Publication No. 85-23, revised 1996) and with UCLA Policy 990 on the Use of Laboratory Animal Subjects in Research (revised 2010).

Senior Editor

  1. Richard Aldrich, The University of Texas at Austin, United States

Reviewing Editor

  1. José D Faraldo-Gómez, National Heart, Lung and Blood Institute, National Institutes of Health, United States

Reviewers

  1. Blanca Rodriguez, Department of Computer Science, University of Oxford, Oxford, United Kingdom
  2. Colleen Clancy

Publication history

  1. Received: March 16, 2018
  2. Accepted: September 24, 2018
  3. Accepted Manuscript published: September 25, 2018 (version 1)
  4. Version of Record published: October 29, 2018 (version 2)

Copyright

© 2018, Rees et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 478
    Page views
  • 91
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Computational and Systems Biology
    2. Developmental Biology
    Inna Averbukh et al.
    Research Article
    1. Computational and Systems Biology
    2. Neuroscience
    Dilawar Singh, Upinder Singh Bhalla
    Research Article Updated