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

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

Similarly to Shiferaw et al. (2003).

(20) Jup=vupci2ci2+cup2

Na+dynamics

Request a detailed protocol

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

Request a detailed protocol

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)

Request a detailed protocol

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)

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

The sarcolemmal Ca2+ ATPase (IPMCA)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

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

Inward rectifier K+ current (IK1)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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)

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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

Request a detailed protocol

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).

Data availability

Gene expression data has been deposited in GEO under accession code GSE48760

The following previously published data sets were used
    1. Rau CD
    2. Wang J
    3. Wang Y
    4. Lusis AJ
    (2013) NCBI Gene Expression Omnibus
    ID GSE48760. Transcriptomes of the hybrid mouse diversity panel subjected to Isoproterenol challenge.

References

  1. Book
    1. Bers D
    (2001)
    Excitation-Contraction Coupling and Cardiac Contractile Force, 237
    Springer Science & Business Media.
    1. Delbridge L
    2. Satoh H
    3. Yuan W
    4. Bassani J
    5. Qi M
    6. Ginsburg KS
    7. Samarel AM
    8. Bers DM
    (1997)
    Cardiac myocyte volume, Ca2+ fluxes, and sarcoplasmic reticulum loading in pressure-overload hypertrophy
    American Journal of Physiology-Heart and Circulatory Physiology 272:H2425–H2435.
    1. Fox JJ
    2. McHarg JL
    3. Gilmour RF
    (2002) Ionic mechanism of electrical alternans
    American Journal of Physiology-Heart and Circulatory Physiology 282:H516–H530.
    https://doi.org/10.1152/ajpheart.00612.2001

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).

Version 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

  • 2,359
    Page views
  • 208
    Downloads
  • 18
    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)

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

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

  1. Colin M Rees
  2. Jun-Hai Yang
  3. Marc Santolini
  4. Aldons J Lusis
  5. James N Weiss
  6. Alain Karma
(2018)
The Ca2+ transient as a feedback sensor controlling cardiomyocyte ionic conductances in mouse populations
eLife 7:e36717.
https://doi.org/10.7554/eLife.36717

Share this article

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

Further reading

    1. Cancer Biology
    2. Computational and Systems Biology
    Bingrui Li, Fernanda G Kugeratski, Raghu Kalluri
    Research Article

    Non-invasive early cancer diagnosis remains challenging due to the low sensitivity and specificity of current diagnostic approaches. Exosomes are membrane-bound nanovesicles secreted by all cells that contain DNA, RNA, and proteins that are representative of the parent cells. This property, along with the abundance of exosomes in biological fluids makes them compelling candidates as biomarkers. However, a rapid and flexible exosome-based diagnostic method to distinguish human cancers across cancer types in diverse biological fluids is yet to be defined. Here, we describe a novel machine learning-based computational method to distinguish cancers using a panel of proteins associated with exosomes. Employing datasets of exosome proteins from human cell lines, tissue, plasma, serum, and urine samples from a variety of cancers, we identify Clathrin Heavy Chain (CLTC), Ezrin, (EZR), Talin-1 (TLN1), Adenylyl cyclase-associated protein 1 (CAP1), and Moesin (MSN) as highly abundant universal biomarkers for exosomes and define three panels of pan-cancer exosome proteins that distinguish cancer exosomes from other exosomes and aid in classifying cancer subtypes employing random forest models. All the models using proteins from plasma, serum, or urine-derived exosomes yield AUROC scores higher than 0.91 and demonstrate superior performance compared to Support Vector Machine, K Nearest Neighbor Classifier and Gaussian Naive Bayes. This study provides a reliable protein biomarker signature associated with cancer exosomes with scalable machine learning capability for a sensitive and specific non-invasive method of cancer diagnosis.

    1. Computational and Systems Biology
    2. Immunology and Inflammation
    Alain Pulfer, Diego Ulisse Pizzagalli ... Santiago Fernandez Gonzalez
    Tools and Resources

    Intravital microscopy has revolutionized live-cell imaging by allowing the study of spatial–temporal cell dynamics in living animals. However, the complexity of the data generated by this technology has limited the development of effective computational tools to identify and quantify cell processes. Amongst them, apoptosis is a crucial form of regulated cell death involved in tissue homeostasis and host defense. Live-cell imaging enabled the study of apoptosis at the cellular level, enhancing our understanding of its spatial–temporal regulation. However, at present, no computational method can deliver robust detection of apoptosis in microscopy timelapses. To overcome this limitation, we developed ADeS, a deep learning-based apoptosis detection system that employs the principle of activity recognition. We trained ADeS on extensive datasets containing more than 10,000 apoptotic instances collected both in vitro and in vivo, achieving a classification accuracy above 98% and outperforming state-of-the-art solutions. ADeS is the first method capable of detecting the location and duration of multiple apoptotic events in full microscopy timelapses, surpassing human performance in the same task. We demonstrated the effectiveness and robustness of ADeS across various imaging modalities, cell types, and staining techniques. Finally, we employed ADeS to quantify cell survival in vitro and tissue damage in mice, demonstrating its potential application in toxicity assays, treatment evaluation, and inflammatory dynamics. Our findings suggest that ADeS is a valuable tool for the accurate detection and quantification of apoptosis in live-cell imaging and, in particular, intravital microscopy data, providing insights into the complex spatial–temporal regulation of this process.