Mathematical modeling of Ca2+ dynamics in the heart has the potential to provide an integrated understanding of Ca2+-handling mechanisms. However, many previous published models used heterogeneous experimental data sources from a variety of animals and temperatures to characterize model parameters and motivate model equations. This methodology limits the direct comparison of these models with any particular experimental data set. To directly address this issue, in this study, we present a biophysically based model of Ca2+ dynamics directly fitted to experimental data collected in left ventricular myocytes isolated from the C57BL/6 mouse, the most commonly used genetic background for genetically modified mice in studies of heart diseases. This Ca2+ dynamics model was then integrated into an existing mouse cardiac electrophysiology model, which was reparameterized using experimental data recorded at consistent and physiological temperatures. The model was validated against the experimentally observed frequency response of Ca2+ dynamics, action potential shape, dependence of action potential duration on cycle length, and electrical restitution. Using this framework, the implications of cardiac Na+/Ca2+ exchanger (NCX) overexpression in transgenic mice were investigated. These simulations showed that heterozygous overexpression of the canine cardiac NCX increases intracellular Ca2+ concentration transient magnitude and sarcoplasmic reticulum Ca2+ loading, in agreement with experimental observations, whereas acute overexpression of the murine cardiac NCX results in a significant loss of Ca2+ from the cell and, hence, depressed sarcoplasmic reticulum Ca2+ load and intracellular Ca2+ concentration transient magnitude. From this analysis, we conclude that these differences are primarily due to the presence of allosteric regulation in the canine cardiac NCX, which has not been observed experimentally in the wild-type mouse heart.
- calcium transient
- C57BL/6 mouse
- Na+/Ca2+ exchanger
the rhythmic contraction of the heart is dependent on tightly regulated intracellular events triggered in a concerted fashion by electrical stimulation. Among a number of regulatory steps, Ca2+ plays, arguably, the most fundamental role in orchestrating the excitation-contraction coupling process. The entry of Ca2+ via voltage-dependent opening of Ca2+ channels during the depolarization phase of the cardiac action potential (AP) triggers the release of Ca2+ [Ca2+-induced Ca2+ release (CICR)] from the sarcoplasmic reticulum (SR) and leads to a global increase in the intracellular Ca2+ concentration ([Ca2+]i), initiating mechanical contractions. Relaxation occurs as [Ca2+]i decays to resting levels, via Ca2+ uptake into the SR through sarco(endo)plasmic reticulum Ca2+-ATPase (SERCA) and Ca2+ extrusion through the Na+/Ca2+ exchanger (NCX) and sarcolemmal Ca2+-ATPase [plasma membrane Ca2+-ATPase (PMCA)].
Due to recent advances in recombinant DNA technology, genetically modified mice have become a popular tool for understanding the mechanisms leading to cardiac diseases, especially those associated with dysfunctional Ca2+ handling (34, 59). With the large amount of experimental data in the mouse now becoming available, mathematical modeling has the potential power to rationalize the data and provide an integrated understanding of the mouse heart. Previously, Bondarenko and colleagues (6) developed a biophysically based electrophysiology model of murine left ventricular (LV) myocytes. However, several limitations with this model (see below) preclude its direct use to quantitatively study Ca2+ handling in the transgenic (TG) mouse heart. Therefore, the focus of this study was to develop a species- and temperature-consistent model of Ca2+ dynamics in LV myocytes of the C57BL/6 mice, with the aim of quantitatively capturing the experimentally observed Ca2+-handling properties. Simulation results were validated against the experimentally observed frequency response of Ca2+ dynamics, AP shape, dependence of AP duration (APD) on cycle length, and electrical restitution.
This model was then applied to analyze the role of NCX in cardiac ventricular myocytes, which acts as the main mechanism for sarcolemmal Ca2+ extrusion and competes with SERCA to lower [Ca2+]i during relaxation. Altered NCX activity and expression have been associated with heart failure, such as in the case of cardiac hypertrophy [for a review, see Reuter et al. (45)]. Using TG mouse models, the effects of overexpression of NCX have been studied by several groups (1, 44, 52, 53, 59). In ventricular myocytes from mice with heterozygous overexpression of the canine cardiac NCX, an unchanged (1, 52, 53) or a slightly higher (59) [Ca2+]i transient magnitude, a faster rate of decay of the [Ca2+]i transient, and an unchanged (1, 59) or enhanced (52, 53) SR Ca2+ content have often been reported. Furthermore, allosteric regulation of NCX by [Ca2+]i, which appears to be absent in the wild-type (WT) mouse heart, has been observed in these TG mice (56). We hypothesized that the mechanisms underlying the observed changes in global Ca2+ handling and the significance of allosteric regulation in the TG murine heart can be better understood through an integrated approach using mathematical modeling. Thus, the second aim of this study was to use the model developed to study the implications of NCX overexpression in cardiac function by simulating and comparing Ca2+ dynamics in the NCX-overexpressing murine heart with and without allosteric regulation.
Following the approach of Niederer et al. (41), to investigate the experimental data dependencies of the Bondarenko et al. model (6), we found significant heterogeneities in the data sources used to determine its model parameters (see Supplemental Material, appendix a, Fig. A1).1 While many of the ion channel formulations, including the fast Na+ current (INa) and various K+ currents [rapid inactivating transient outward K+ current (IKto,f), slow inactivating transient outward K+ current (IKto,s), rapidly activating slowly inactivating K+ current (IKur), and noninactivating steady-state K+ current (IKss), etc.] were derived from murine LV myocytes, the intracellular Ca2+-handling system was based mostly on the Jafri et al. (26) model for guinea pigs. No direct experimental measurements in the mouse have been used to refit these parameters, and very little explanation has been given as to how the parameter values were chosen. Furthermore, the Bondorenko et al. (6) model characterizes the murine AP at 25°C. However, to simulate murine cardiac myocytes at more physiological conditions requires the reparameterization of the model to capture the faster ion channel kinetics at higher temperatures. Recently, more data have become available that now allows a model at a physiologically relevant temperature to be developed.
The limitations of the Bondarenko et al. model, as explained above concerning its heterogeneous sources of experimental data and the temperature difference, were the main motivations for the development of the present model. In additional, some properties of simulated Ca2+ dynamics and AP kinetics were found to be inconsistent with experimental observations, as explained below.
It is also worth noting that the original formulation of the Bondarenko et al. model does not account for the ionic concentration changes introduced by the stimulus current, which is required for charge conservation and to avoid gradual drift in the variables (11, 24). Following the approach of Hund et al. (24), we assigned the stimulus current to be a K+ carrier in all the simulations discussed below.
The time to half decay of the simulated [Ca2+]i transients using the Bondarenko et al. model paced at a lower pacing frequency (0.5 Hz, 49 ms) is <50 ms, which is significantly shorter than that reported experimentally (175 ± 25 ms) at the same pacing frequency and under similar extracellular ion concentrations (17). In addition, at 1 Hz, simulated Ca2+ flux through PMCA is similar in magnitude to that through SERCA and significantly greater than that through NCX. The percent contribution of SERCA, NCX, and PMCA to total Ca2+ removal is 46.1%, 6.5%, and 47.4%, respectively, which is not consistent with experimental observations in the mouse showing that SERCA is responsible for >90% of total Ca2+ removal and PMCA is responsible for <1% (34). Finally, simulated AP using the Bondarenko et al. model does not reproduce APD restitution properties, i.e., decreasing APD with decreasing diastolic interval. Such a phenomenon has been observed experimentally in the intact mouse heart (30) and is believed to be an important factor underlying cardiac arrhythmia (18).
Isolation of Myocytes
LV myocytes were isolated using a standard enzymatic dispersion technique (46). Briefly, the heart was initially perfused for 3 min with a nominally Ca2+-free solution followed by a further 8–9 min with an enzyme-containing solution (collagenase, 1 mg/ml, Worthington Biochemical, Ca2+: 0.05 mmol/l). Unless specified otherwise, myocytes were superfused with a modified Tyrode solution containing (in mmol/l) 134 NaCl, 5.4 KCl, 1.2 MgCl2, 1.4 CaCl2, 11.1 glucose, and 5 HEPES at pH equal to 7.4 (adjusted with NaOH). The LV free wall was isolated and placed in a separate flask containing fresh enzyme solution. Myocytes were harvested after further 5- and 10-min digestion periods, washed twice, and resuspended in storage solution containing (in mmol/l) 120 NaCl, 5.4 KCl, 5 MgSO4, 0.2 CaCl2, 5 Na-pyruvate, 20 glucose, 20 taurine, and 10 HEPES at pH equal to 7.4 (adjusted with NaOH). The myocyte suspension was stored at room temperature, and cells were used within 8 h of isolation. All protocols were in accordance with Home Office (UK) guidelines and have been approved by Oxford University's Ethics Committee in accordance with Schedule I Procedures.
Measurement of [Ca2+]i Transients
The [Ca2+]i transient was evaluated in fura-2 AM-loaded (5 μmol/l, Molecular Probes) LV myocytes field stimulated at 0.5, 1, and 3 Hz or exposed to caffeine (10 mmol/l for 10 s) after the SR was loaded by a train of >50 stimuli at 3 Hz. Although the basal frequency of 3 Hz is slower than the typical heart rate of a mouse (10 Hz), it is the highest frequency at which synchronized [Ca2+]i2+]i transient was calculated as the difference between diastolic and peak Ca2+ fluorescence. [Ca2+]i transients from individual myocytes were obtained by averaging the recordings over a stable period of ∼10 s to provide characteristic [Ca2+]i transients for analysis.
Measurement of Ca2+ Current
L-type Ca2+ current (ICaL) was measured using the whole cell configuration of the patch-clamp technique (Axopatch 200A, Axon Instruments) in the presence of 1 mmol/l 4-aminopyridine (4-AP) to block fast transient outward K+ currents. ICaL was elicited by 200-ms step depolarizations from a holding potential of −40 mV to a range of test potentials at a stimulation frequency of 1 Hz. Peak inward current was measured with respect to the current at −40 mV. The average cell membrane capacitance, measured by applying a −10-mV pulse from a holding potential of −40 mV for 50 ms, was 138.66 pF. All experiments were conducted at 35°C.
The ratiometric fluorescence was converted into free [Ca2+]i following the method of Chorvatova and Hussain (10) using the following equation: (1) where Kd is the dissociation constant (0.251 μM), R is the fluorescence ratio recorded at the two excitation wavelengths (fluorescence at 365 and 380 nm), Rmin (0.909) and Rmax (2.695) are the fluorescence ratios under Ca2+-free and Ca2+-saturating conditions, respectively, and Fmax380 (326.286) and Fmin380 (83.543) are the maximum and minimum fluorescence intensities excited at 380 nm.
MATLAB (The MathWorks, Matick, MA) was used for all data analysis and parameter fitting. All error bars shown represent SEs. For parameter fitting, we used either a subspace trust region method (“lsqcurvefit” function) or the unconstrained nonlinear optimization function based on the Nelder-Mead simplex method (“fminsearch” function).
To facilitate transparent and unambiguous dissemination of the results of this study, we have made the model available on the CellML website (www.cellml.org), along with metadata links to the experimental data sets used for parameterization of the Ca2+-handling components. A full explanation of the data sets and how they are linked with each model component is also provided.
Ca2+ flux through NCX (WT murine).
Current through the NCX (INCX) was modeled using the following formulation from Weber et al. (56): (2) where V is the transmembrane voltage, VNCXmax is maximum voltage through NCX, [Na+]i is intracellular Na+ concentration, [Ca2+]o is extracellular Ca2+ concentration, η is the energy barrier position controlling the voltage dependence of INCX, F is Faraday's constant (96.5 C/mmol), R is the universal gas constant (8.314 J·mol−1·K−1), and T is the absolute temperature (308 K). KmAllo is the affinity for allosteric [Ca2+]i; KmCao, KmCai, KmNao, and KmNai are the transport affinities for Ca2+ and Na+ outside and inside. ksat (0.27) assures saturation at very negative voltage. KmAllo was set to zero as allosteric regulation of the exchanger by [Ca2+]i was not observed in WT mouse ventricular myocytes (56). The relationship between INCX and the flux of Ca2+ through the NCX (JNCX) into the cytosol was defined as follows: (3) where Cm is the average cell membrane capacitance (138.66 pF), and the myoplasmic volume of the murine ventricular myocyte (Vmyo) was set at 22 × 10−6 μl (42).
Caffeine-induced SR Ca2+ release is a commonly used experimental protocol to investigate the function of NCX. The time courses of caffeine-induced [Ca2+]i transients from individual myocytes (n = 18) were first averaged and converted to free [Ca2+]i using Eq. 1. [Ca2+]i was then converted to the total cytosolic Ca2+ concentration ([Ca2+]total) using the following fast buffering approximation proposed by Trafford and Eisner (54): (4) where Kd,buffer (0.60 μM) is the affinity of the buffer, Bmax (109 μM) is the concentration of the buffer, and a is an arbitrary offset (54). Since this equation represents the total amount of cytosolic Ca2+ buffering, the original Ca2+-buffering equations in the Bondarenko et al. model (Eqs. A16 and A17 in Ref. 6) were removed in our new formulation. The time derivative of [Ca2+]total in Eq. 4 determines the net flux of Ca2+ across the sarcolemma (Jsarcolemma). Assuming that the contribution of PMCA to the decay is comparatively small (34), we attributed 90% of the total Ca2+ extrusion to NCX and 10% to PMCA (34). The flux of Ca2+ through PMCA (JPMCA) takes the same form as that used by Bondarenko et al. (6). The fitted values of maximum pump rate (IPMCAmax) and affinity for Ca2+ (Km,PMCA) are 0.096 pA/pF and 0.289 μM, respectively. The [Ca2+]i dependence of JNCX was fitted according to Eqs. 2 and 3 as follows.
To determine [Na+]i, parameters for the Na+/K+ pump were first adjusted according to experimental measurements of the pump current at 35°C by Berry and colleagues (4). Specifically, the Hill coefficient (nH) was set to 2.4, and the affinity of the pump to Na+ was adjusted to 16.6 mM. IPMCAmax was then determined to be 2.486 pA/pF, such that Na+/K+ pump current amplitude at 0 mV in 5 mM intracellular K+ concentration and 50 mM [Na+]o was 1.4 pA/pF. These new parameter values were incorporated into the AP framework paced at 1 and 3 Hz to obtain [Na+]i values of 10 and 12.5 mM, respectively. [Na+]i in quiescent murine myocytes has been reported to be between 14.9 and 17.1 mM at 25°C (59). However, the activity of the pump is expected to be enhanced at higher temperatures, resulting in a lower [Na+]i, meaning that our estimates remain consistent with cellular functioning at 35°C. In addition, we also assumed that the transmembrane potential remained constant during caffeine application at −80 mV. The validity of this assumption was confirmed as the experimental protocol was reproduced in our model, i.e., the simulated transmembrane potential was −79.1 mV during caffeine application. For parameterization, [Na+]i was set at 12 mM. Extracellular ion concentrations, set to be consistent with the experimental setup, were 134 mM [Na+]o and 1.4 mM [Ca2+]o. Other constants included the following: F (96.5 C/mmol), R (8.314 J·mol−1·K−1), T (308 K), Cm (138.66 pF), and Vmyo (22 × 10−6 μl).
Initially for the fit, only VNCXmax was allowed to vary, and we set all other parameters to the same values as those used by Weber et al. (56). Figure 1A, top, shows a caffeine-induced [Ca2+]i transient, expressed as the fura-2 AM ratio, from a single myocyte with prior pacing at 3 Hz. The average time course of the transient, expressed as [Ca2+]i, is shown in Fig. 1A, bottom. [Ca2+]i was then converted to [Ca2+]total using Eq. 4. Assuming that no Ca2+ is taken up by the SR and that the SR is depleted during the decay phase of the caffeine transient, the time derivative of [Ca2+]total is equal to Jsarcolemma. After taking into account the small background Ca2+ leak current and assuming that the contribution of PMCA to the decay is comparatively small (34), we attributed 90% of the total Ca2+ extrusion to NCX. The dependence of JNCX on [Ca2+]i can then be shown by plotting JNCX against the corresponding [Ca2+]i at each point in time, as shown in Fig. 1C. The simulated JNCX as a function of [Ca2+]i, with the exchanger rate reparameterized (VNCXmax: 3.939 pA/pF), is shown as the solid line in Fig. 1C. The root mean square deviation (RMSD) between the predicted JNCX and the values calculated from experimental data was 1.219 × 10−3.
In the original formulation proposed by Weber et al. (56), there are a total of seven parameters. Taking into account the Haldane relationship (56), i.e., KmCai·KmNao3 = KmNai3·KmCao and hence KmNao3 = KmNai3·KmCao/KmCai, there are now six free parameters, as shown in Table 1. We then carried out a sensitivity analysis on these parameters following the approach by Fink et al. (16). Briefly, we defined variable y as the RMSD between the predicted and experimentally calculated JNCX, representing the quality of the fit, and y = y(p), where p is the parameter that we are interested in. It follows that at p = p0, y0 = y(p0), and at p = p0 + Δp, we define Δy = y(p0 + Δp) − y(p0). The sensitivity (σy,p) of y with respect to parameter p at p0 was defined as follows: (5) By setting p0 to be the original parameter value used by Weber et al. (except for VNCXmax, which was set to be the fitted value) and Δp to be 5% of p0, the sensitivity of RMSD to each parameter (listed in Table 1) were calculated using this method.
It can be seen that the quality of the fit is most sensitive to VNCXmax and moderately sensitive to KmCai and ksat. It is insensitive to η, KmCao, and KmNai under this protocol. To check whether allowing more parameters to vary would give a better fit, we tried fitting the three most sensitive parameters, KmCai, ksat, and VNCXmax, yielding an RMSD of 1.211 × 10−3. The quality of fit was only marginally better than that obtained by fitting VNCXmax only (RMSD: 1.219 × 10−3), and the differences between the fits were not visually detectable (not shown). For this reason, we did not vary the values for KmCai and ksat in the final model. Thus, Eq. 2 can be simplified by setting parameters KmCao and KmNai to be constants and making use of the following relation: KmNao3 = KmNai3·KmCao/KmCai. Parameter η was kept as a variable as it governs the voltage dependence of INCX. The simplified formulation can be written as follows: (6)
Ca2+ flux through SERCA.
The [Ca2+]i dependence of Ca2+ flux through SERCA (JSERCA) was determined by integrating the NCX model fitted above with the decay phase of the [Ca2+]i transients elicited by field stimulation at a range of frequencies. During a field-stimulated [Ca2+]i transient, the sources of [Ca2+]i are assumed to be ICaL and SR Ca2+ release, whereas the [Ca2+]i sinks are SERCA, NCX, and PMCA. The inactivation time constant for the L-type Ca2+ channel is faster than 40 ms (see below), whereas CICR was assumed to have a duration of no more than 80 ms based on a study (49) performed in rats. Therefore, the effects of Ca2+ influx through ICaL and Ca2+ release through ryanodine receptors (RyRs) can be neglected when analyzing the later part of the [Ca2+]i transient. To characterize SERCA, Ca2+ fluxes were calculated after 70 ms (3 Hz) or 90 ms (0.5 and 1 Hz) had elapsed from the peak of the [Ca2+]i transient.
At each frequency, the ensemble flux of Ca2+ (Jtotal) during the decay of the transient was calculated using the same method as previously described for the caffeine-induced [Ca2+]i transient. Jsarcolemma, estimated from the NCX and PMCA models derived above, was then subtracted from Jtotal to give the net flux of Ca2+ into the SR (JSERCA). Here, we relied on the assumption that the Ca2+ efflux due to mitochondria uptake is minimal under physiological conditions (5). To summarize, JSERCA was obtained from the following: (7) In the present study, JSERCA was assumed to be made up of two components: SERCA uptake (Jup) and leak (Jleak ). The formulation of the Jup model was adopted from Jafri et al. (26), giving the following equations: (8) (9) where the negative sign denotes Ca2+ movement from the intracellular space to the SR, Vup is the maximum pump rate of SERCA at saturating [Ca2+]i, Km,up is the Ca2+ affinity of SERCA, and Hill coefficient nH (equal to 2) was set to be consistent with experimental findings (36).
Jleak was, in turn, defined as follows: (10) where the SR leak rate (Vleak) determines the dependence of Jleak on the Ca2+ concentration gradient between the network SR and the cytosol [Ca2+] concentration in the network SR ([Ca2+]NSR) − [Ca2+]i.
Figure 2A, top, shows an example of experimental recordings from a myocyte field stimulated at 3 Hz. The average [Ca2+]i transient from individual myocytes (n = 22) is shown in Fig. 2A, bottom. Jtotal during the decay of the [Ca2+]i transient was calculated using the same method as previously described for the caffeine-induced [Ca2+]i transient. Jsarcolemma, estimated from the NCX and PMCA models derived above, was then subtracted from Jtotal to give JSERCA. JSERCA was then used to reparameterize Km,up and Vup and provide an estimate of Jleak. Figure 2B shows the experimentally derived total Ca2+ flux (points) compared with the fitted total flux (dotted line) consisting of JNCX, JPMCA, Jup, and Jleak. Vleak was then determined in a way such that in whole cell simulations at 3 Hz, Jleak was close to that estimated from parameter fitting. Jleak, Vup, and Km,up were fitted to JSERCA at 3 Hz, yielding Jleak equal to 0.037 μmol·l cytosol−1·ms−1, Vup equal to 0.82 μM/ms, and Km,up equal to 0.412 μM.
An important assumption underlying the Jafri et al. (26) formulation of SERCA is that Vup and Km,up remain constant across a range of pacing frequencies. In our model, fitting Jup to experimental data recorded at different pacing frequencies suggested that Jup is upregulated at 3 Hz compared with 0.5 and 1 Hz (n = 22). Figure 2C, top, shows an indicative example of an experimental recording of the frequency response of [Ca2+]i transients from one myocyte. Figure 2D shows the frequency-dependent changes in the Ca2+ dependence of Jup, with Jup reaching a greater magnitude (more negative) for a given level of [Ca2+]i at 3 Hz compared with lower frequencies. The values for Vup at 0.5 and 1 Hz were determined by fitting Jup according to Eq. 9. This yielded Vup values of 0.55 and 0.65 μM/ms for 0.5 and 1 Hz, respectively. Existing literature supports the finding that the Ca2+ uptake rate through SERCA becomes faster at higher frequencies to facilitate frequency-dependent acceleration of relaxation (23). Frequency-dependent acceleration of relaxation is thought to involve Ca2+/calmodulin kinase II (CaMKII) through the phosphorylation of key Ca2+ regulatory proteins such as phospholamban (32) and SERCA (57), although the exact mechanism remains debated (21, 34). Picht and colleagues (43) recently carried out a quantitative analysis of cytosolic Ca2+ removal during [Ca2+]i decline and found that in mice, increases in the level of active CaMKII at higher pacing frequencies lead to the elevation of SERCA Ca2+ uptake by way of accelerating Vup without affecting its affinity.
To model the frequency dependence of Jup, the CaMKII regulatory pathway was incorporated into our model using the formulation proposed by Hund and Rudy (HRd) for canine ventricular myocytes (25). Many of the parameter values used in the HRd model were taken from a mathematical model of CaMKII function in the brain, as experimental data for its neuronal isoform are much more abundant than that for the cardiac isoform. We decided to use almost identical parameter values as those in the HRd model, since there are no sufficient data available, to our knowledge, to allow reparameterization. While highly simplified in its description of CaMKII dynamics and thus computationally efficient, the HRd model is able to produce steady-state [Ca2+]i transients for a wide range of pacing frequencies that agree with experimental measurements (25).
The state variables fa, fb, and ft in our model correspond to CaMKactive, CaMKbound, and CaMKtrap in the HRd model, respectively. Parameters f0, α, and β correspond to CaMK0, αCaMK, and βCaMK in the HRd model, respectively. Specifically, the fractional occupancies of bound CaMKII (fb) and trapped CaMKII (ft) are described below, and the total fractional occupancy of active CaMKII (fa) is given by the sum of fb and ft, as follows (11) (12) (13) where f0 (5%) is the fractional occupancy at equilibrium of the low-affinity calmodulin-binding sites and Km,CaM (0.7 μM) is the concentration of [Ca2+]i for half-maximal activation of the kinase (60). α (0.05 ms−2) is the autophosphorylation rate, and β, assumed to be 2.0 × 10−4 ms−1, is the rate of dephosphorylation. Minor modifications to the HRd model for CaMKII include, first, that fb is a function of [Ca2+]i instead of the steady-state Ca2+ concentration. This was done because we wanted to use experimentally measured [Ca2+]i transients as an input to simulate the level of active CaMKII. Second, the value for Km,CaM was changed from the original value of 1.5 to 0.7 μM, as a more recent study (60) suggested half-maximal activation of CaMKII occurs at 0.5–1 μM free Ca2+. Finally, the value for β was changed from 0.00068 to 0.0002 ms−1 to avoid large oscillations in fa over each cardiac cycle, which was seen when the original value was used in our simulations.
The dependence of Jup on CaMKII was then modeled by defining Vup from Eq. 9 as a function of fa, as shown in Fig. 2C, bottom: (14) (15) (16) where V̄up is the maximal uptake rate in the absence of CaMKII, denotes the maximum CaMKII-dependent increase in Vup, Km,CaMK is the half-saturation coefficient and n is Hill coefficient. The line of best fit for Vup plotted over fa for each frequency yielded V̄up (0.396 μM/ms), (2.998), Km,CaMK (1.244) and n (2.583).
An example of the experimentally recorded transmembrane currents elicited by test potentials between −30 and 20 mV is shown in Fig. 3A. INa was inactivated by applying a holding potential of −40 mV, and the transient outward K+ currents were blocked due to the presence of 4-AP (see Measurement of Ca2+ Current). Therefore, the peak transmembrane current measured at each test potential, as shown in Fig. 3B, can be taken as a good approximation of ICaL. This assumption was validated using the final model (as shown in the Supplemental Material, appendix a).
Under the particular experimental protocol used, the decay phase of the transmembrane current consists of both ICaL and the late K+ currents, such as the slow delayed rectifier K+ current (IKs) and rapid delayed rectifier K+ current (IKr), which are insensitive to 4-AP. It was therefore not possible to directly acquire the inactivation kinetics of the channel. As an approximation, we used the formulations proposed by Bondarenko et al., with reparameterized parameter values (see below), to approximate the magnitudes and time courses of the late K+ currents and subtracted them from the experimentally recorded transmembrane current. The remaining current was assumed to consist of ICaL only and was fitted to a single exponential function to give the time constant of decay (τ) at each test potential, shown in Fig. 3C.
The original formulation for the L-type Ca2+ channel in the Bondarenko et al. (6) model consists of eight channel states with multiple pathways between different states. The complexity of this formulation is problematic given that only a limited set of experimental data are available to parameterize the model. We therefore used a simpler formulation, with a smaller number of parameters, which was still able to reproduce our experimental observations.
We chose a modified version of the L-type Ca2+ channel model developed by Hinch et al. (22) for rat ventricular myocytes, which is a simplification of the Jafri et al. (26) model. In this formulation, transitions from the closed state C to the open state O and vice versa are governed by voltage-dependent transition rates α+ and α−, respectively. Transition from the closed state C to the inactivated state I is a function of Ca2+ concentration in the dyadic space ([Ca2+]ds), which reflects its Ca2+-dependent nature. The voltage-dependent inactivation gate y in the original Jafri et al. (26) model was retained here to capture the voltage inactivation kinetics of the channel. The transition rates were defined the same way as proposed by Hinch et al. (22) with nine parameters (see Supplemental Material). The following formulation for the voltage-dependent inactivation gate y was similar to that proposed by Jafri et al. (26): (17) which has nine parameters (C1–C9), yielding a total of 18 parameters. The flux of Ca2+ through the L-type Ca2+ channel (JCaL) into the dyadic space is then governed as follows: (18) where P̄CaL is the permeability of the channel, O is the channel open probability, and δ = zF/RT, with z being the valance of Ca2+, F being 96.5 C/mmol, R being 8.314 J·mol−1·K−1. This yields the following equation for ICaL: (19) where Vds is the volume of the dyadic space, which was set to 2.2 × 10−8 μl to keep the same myoplasmic and dyadic space volume ratio as that in the Bondarenko et al. model, and Acap is capacitive membrane area.
Hinch et al. (22) fitted the potential when half of the L-type Ca2+ channels are open (VL) as well as the width of opening potentials (ΔVL) to the current-voltage (I-V) curve of the current (22). However, it is likely that the other 16 parameters also need to be reassessed due to differences between the rat and mouse. Therefore, the parameter values were identified as follows.
First, some of the parameters governing the voltage-dependent inactivation can be determined from experimental data obtained with Ba2+ as the charge carrier for the L-type Ca2+ channel (40). Using double-pulse voltage-clamp protocols (as shown in Supplemental Material, appendix a, Fig. A4), Nguemo et al. (40) measured the availability-voltage relationship of the current in the presence of Ba2+ in adult mouse ventricular myocytes at 37°C. The above measurements were used to fit parameters C1, C2, and C7, yielding C1 = 33 mV, C2 = 8.23 mV, and C3 = 315 ms.
Second, we carried out a sensitivity analysis, similar to that used for the NCX model, to examine the importance of the other 13 parameters in determining the quality of the fit to the experimentally measured I-V curve and the experimentally calculated voltage-dependent inactivation time constant (τ) of the current. As shown in Table 2 (see also Supplemental Material, appendix a, Fig. A4), the I-V curve and the inactivation kinetics were both very sensitive to the following parameters: proportion of time closed in the open mode (ΦL), biasing to make inactivation a function of voltage (a), PCaL, and time switching between closed and open states (tL). Further visual investigation (not shown) revealed that the parameters ΦL and tL determine the slope of the I-V curve, whereas the parameters a and PCaL only scale the I-V curve. The inactivation kinetics are also very sensitive to concentration at inactivation (KL), inactivation time constant (τL), C6, C8, and C9. Both the I-V curve and the inactivation kinetics are insensitive to the following parameters: biasing to make inactivation a function of voltage (b), C3, C4, and C5.
Therefore, we fitted PCaL (2.5 ms−1), ΦL (1.798), and tL (1.5), along with VL (0 mV) and ΔVL (6.449 mV), to the I-V curve, as shown in Fig. 3B, with the inclusion of the NCX and SERCA models fitted previously and the conductivities of the K+ currents set to zero. Parameters KL (0.3 μM), τL (1,150 ms), C6 (10 ms), C8 (30 mV), and C9 (4 mV) were then fitted to the experimentally calculated voltage-dependent inactivation time constant (τ) of the current, as shown in Fig. 3C. In addition, parameter b (0.4) was chosen such that the magnitude of the current elicited by field stimulation at 3 Hz was ∼50% of that at 0.5 Hz according to experimental observations in mouse ventricular myocytes (2). Finally, parameters C3, C4, and C5 were kept to their original values as those proposed by Jafri et al. (26). Simulated time courses of the current elicited by test potentials between −20 and 20 mV in 10-mV intervals are shown in Fig. 3D.
Other transmembrane currents.
To ensure temperature consistency of the model, formulations of other transmembrane currents in the Bondarenko et al. model, including the various K+ currents and INa, were reparameterized to experimental measurements at a more physiological temperature. There are seven types of K+ currents with different kinetics in the Bondarenko et al. model, six of which are considered here: IKto,f, IKur, IKss, IKr, IKs, and time-independent K+ current (IK1). We chose to exclude IKto,s since it is only found in the septum region of the heart (58) and the cardiomyocytes used in this study (see Isolation of Myocytes) were from the LV free wall. Parameters for IKs were kept the same as those used in the Bondarenko et al. (6) model as its expression in the murine heart was questionable and the magnitude was relatively small. Whenever possible, each of these currents was parameterized individually using the appropriate experimental data from literature as explained below. When direct measurements were not available, parameters were fitted to the ensemble K+ current (14) and validated using experimentally measured APDs.
For IKto,f, the kinetics of the inactivation state (ito,f) can be rewritten as follows: (20) where ito,f,ss is the proportion of the channel that is inactivated at steady state. The inactivation time constant (τito,f) was fitted to the experimentally measured time course of inactivation shown in Fig. 4A, yielding (21) ito,f,ss was fitted to the experimentally measured voltage dependence of steady-state inactivation, as shown in Fig. 4B, yielding (22) The voltage dependence of the activation of IKto,f was fitted to experimentally measured transient and sustained components of the ensemble outward K+ currents (14) as follows: (23) where αa and βa are the voltage-dependent rate constants for activation. For IKur, the proportion of the channel that is inactivated at steady state (iss) was fitted to the experimentally measured voltage dependence of steady-state inactivation corrected for the presence of divalent ions (Co2+), as shown in Fig. 4C. The time constant of inactivation (τiur) was fitted to the experimentally measured recovery from steady-state inactivation, as shown in Fig. 4D. Other parameters were the same as those in the Bondarendo et al. model (6). (24) For IKr, the conductance of the channel and the rate constants from state CK2 to state OK (αa1) and from state OK to state IK (αi) were fitted to the experimentally measured I-V relationship (33), as shown in Fig. 4E: (25) For IK1, the conductance of the channel (0.35 mS/μF) was fitted to the experimentally recorded I-V relationship (29), as shown in Fig. 4F.
Finally, experimentally measured transient and sustained components of the ensemble outward K+ currents (14), as shown in Fig. 5, A and B, were used to parameterize the conductance of IKto,f (0.535 mS/μF), IKur (0.250 mS/μF), and IKss (0.060 mS/μF). These were also used to fit the voltage dependence of variable ass, which governs the steady-state activation of IKur and IKss, as follows: (26) Due to its large magnitude and fast kinetics, experimental data on the activation and inactivation kinetics of INa done at physiological temperatures are very limited, and we did not fit INa explicitly. To determine the scope of variations resulting from uncertainty in INa, a sensitivity analysis of AP shape to the conductance of INa (GNa) was carried out. A 50% increase in GNa caused a 25% increase in the maximum upstroke velocity of the AP (Vmax), a 10% increase in AP overshoot, a 3% decrease in APD at 25% repolarization (APD25), a 5% increase in APD at 90% repolarization (APD90), and a 1.5% increase in [Na+]i. Therefore, GNa was adjusted to match the experimentally measured Vmax (177 V/s) (12).
RyR Ca2+ release and intra-SR buffering parameters.
The rate of SR Ca2+ release through RyRs during CICR was fitted to match the experimentally measured time to peak of the [Ca2+]i transient at 3 Hz (25 ± 1 ms). Parameters determining the rate of changes in PRyR [the first and second constants in Eq. A15 in the Bondarenko et al. (6) model, termed APRyR and BPRyR, respectively, in the Supplemental Material, appendix b] were adjusted from 0.04 and 0.1, respectively, to 0.01 and 2.0, respectively, such that the peak opening probability of RyRs is higher with a shorter opening duration. (27) The total concentration of calsequestrin was set to 50 mM, or 175 μmol Ca2+-binding sites/l cytosol, close to the estimates by Bers (5) of the range between 175 and 350 μmol/l cytosol. The affinity of calsequestrin for Ca2+ was 630 μM according to measurements found in the literature (47) .
Background Ca2+ current.
Background Ca2+ current (ICab) was formulated as linear leakage current as follows: (28) where GCab is the conductance of the channel and ECaN is the Ca2+ reversal potential. GCab (7.0 ×10−4 mS/μF) was set such the simulated diastolic [Ca2+]i level was 0.10 μM at 3 Hz, close to the experimentally measured level (0.11 ± 0.01 μM).
Volumes of intracellular compartments.
With the myoplasmic volume set to be 2.2 × 10−6 μl, the volumes of intracellular compartments were determined by their percentage of the cell volume. Using stereological techniques with electron microscopy, Bossen and colleagues (7) reported the volume fractions of the mitochondria and total SR to be 37.5% and 0.87% of the total cell volume in murine LV myocytes. Assuming the myoplasm to be 62.5% of the total volume, the fractional volume of total SR is thus 1.4% of the myoplasm. In the same study, the volume fraction of the junctional SR was reported to be 0.22% of the total cell volume, equivalent to 25% of the total SR. Accordingly, the volumes of the junctional and network SR were set to be 7.7 × 10−8 and 2.31 × 10−7 μl, respectively. The volume of the dyadic space is reportedly near 0.1% of cell volume (5) and was set to be 2.2 × 10−8 μl.
Canine cardiac NCX overexpression model.
The formulation of NCX in the TG model was based on studies by Weber et al. (56) and Terracciano et al. (52). According to the study of Weber et al. (56), at [Na+]i equal to 10 mM, KmAllo in TG mice was estimated to be 0.15 μM. For this reason, we set KmAllo in our TG model to be 0.15 μM. We then fitted the NCX parameters VNCXmax (15.0 pA/pF) and η (0.6) in the TG model to the I-V relationship measured by Terracciano et al. (52) at room temperature, since to our knowledge no data at 35°C have been published. In their study, INCX was recorded under a voltage-ramp protocol, where the membrane was depolarized from a holding potential of −75 to −40 mV for 200 ms followed by a voltage ramp to 80 mV and back to −120 mV at a velocity of 50 mV/s.
Model Validation (WT Murine)
After adjusting the transmembrane current parameters to a more physiological temperature, the simulated AP had kinetics comparable with experimental observations. At a pacing frequency of 1 Hz (Fig. 5C), the resting potential during diastole was −81 mV, close to the experimental values of −79 and −80 mV (8, 14, 29). The overshoot potential was 38.5 mV compared with 40 and 46 mV, as reported by Brouillette et al. (8) and Knollmann et al. (29), respectively. At 1-Hz pacing frequency, the simulated APD to 50% repolarization (APD50) was 2.92 ms compared with the experimental values of 1.7 ms (27), 2.5 ms (8), 4.7 ms (14), and 6.0 ms (29). APD90 was 19.2 ms, compared with experimental values of 23.6 (20), 25 (29, 35), 15, and 15.4 ms (27). Table 3 shows experimentally measured AP characteristics and the corresponding strains of mice, cell type, and temperature.
Validating the RyR model.
In the proposed model, the properties of the RyR, an important modulator for CICR, were not directly measured. Instead, we adapted the formulation used in the Bondarenko et al. model with some adjustments to parameter values to match the experimentally observed time to peak of the [Ca2+]i transient. Here, we show that our model is able to qualitatively reproduce key characteristics of CICR.
The first set of simulations tested the ability of our model to capture the property of graded release, which refers to the observation that SR Ca2+ release is proportional to the influx of trigger Ca2+. Figure 6A shows the time courses of ICaL (top) and the [Ca2+]i transient (middle) in response to a 200-ms step depolarization from a holding potential of −50 to 0 mV. Figure 6A, bottom, shows the peak ICaL and [Ca2+]i as a function of the test potential between −40 and 40 mV. Peak ICaL showed a voltage dependence with a maximum at 0 mV, and peak [Ca2+]i followed the same qualitative trend.
The second set of model simulations aimed to reproduce the autoregulation observations of Eisner et al. (15), who demonstrated the transient nature of the cell's response to an increase in the open probability of RyRs. In their study, a low concentration of caffeine was applied to field-stimulated ventricular myocytes. This produced a transient increase in the amplitude of the systolic [Ca2+]i transient, which gradually returned to the same steady-state level as that of the control. To validate our model, we repeated the same protocol in our simulations by doubling the open probability of the RyR (PO1 + PO2) during steady-state pacing at 3 Hz. The magnitude of the simulated [Ca2+]i transient, as shown in Fig. 6B, top, increased upon the application of caffeine and gradually decreased over time until it returned to the magnitude before caffeine application. As shown in Fig. 6B, middle, there was a decrease in total SR Ca2+ after caffeine application, which did not recover over time. When examining the net sarcolemmal Ca2+ flux in Fig. 6B, bottom, we found that the loss of Ca2+ from the SR is due to a net Ca2+ efflux immediately after caffeine application resulting from the simultaneous transient increase in the [Ca2+]i transient. As SR Ca2+ decreases, the [Ca2+]i transient also decreases until a new balance point is established, where there is zero net sarcolemmal flux. The above simulation results were consistent with those observed by Eisner et al. (15) and further validate the ability of our model to capture CICR characteristics.
Ca2+ dynamics versus pacing frequency.
As shown in Fig. 7, A and B, simulated [Ca2+]i transients paced at 0.5, 1, and 3 Hz were in good agreement with experimental measurements. Overall, a slightly negative frequency response of [Ca2+]i transients could be observed experimentally, with peak [Ca2+]i decreasing from 0.50 ± 0.09 μM at 0.5 Hz to 0.37 ± 0.06 μM at 3 Hz. This compared well with simulated peak [Ca2+]i values of 0.48 μM at 0.5 Hz and 0.36 μM at 3 Hz. Consistent with our analysis of frequency-dependent SERCA function, the time to half relaxation (RT50) was shortened from 94 ± 6 ms at 0.5 Hz to 65 ± 3 ms at 3 Hz, which closely matched with the simulated values of 96 and 70 ms, respectively. A summary of simulated and experimental [Ca2+]i transients parameters is shown in Table 4.
The time courses of [Ca2+]i fluxes through SERCA, NCX, and PMCA during one cardiac cycle at 3 Hz are shown in Fig. 7C. The percent contribution of SR Ca2+ uptake to total [Ca2+]i removal, calculated as the ratio between the integral of JSERCA over one cardiac cycle and the integral of the sum of the fluxes, was 93.4%. NCX and PMCA contributed 5.9% and 0.7%, respectively.
A simulated caffeine-induced [Ca2+]i transient is shown in Fig. 7D. A train of field stimulations at 3 Hz was applied to the model before caffeine application until the limit cycle had been reached. The effects of caffeine application were replicated by setting the open probability of RyR to 1 and increasing the rate SR Ca2+ leak to 0.1 ms−1, causing a rapid and complete release of Ca2+ from the SR. Experimentally recorded transients from individual myocytes and the average transient are also shown in the same plot for comparison. Peak [Ca2+]i of the simulated transient was 0.49 μM compared with the experimental average of 0.55 μM. The time constant of decay of the simulated transient was 1.27 s compared with 1.36 s recorded experimentally.
Effects of cycle length and restitution.
The effects of pacing frequency on the shape of steady-state APs are shown in Fig. 8, A and B, with pacing frequency increasing from 0.5 to 10 Hz in the direction indicated by the arrow. There were no significant differences in the upstroke and early repolarization phases (APD25) between different frequencies. However, late repolarization, measured by APD at 75% repolarization (APD75) and APD90, was progressively prolonged with increasing cycle length, which matches qualitatively with the experimental characterization of the AP at various cycle lengths in intact murine hearts. It has been suggested that in murine and rat ventricular myocytes, the inward current through the electrogenic NCX during the later phase of the AP is responsible for shaping APD90 (28). To test whether our model agrees with the literature, we eliminated INCX by setting VNCXmax to zero in our model. APs in the absence of INCX over the same range of frequencies as previously indicated are shown in Fig. 8C, with the corresponding APDs in Fig. 8D. It was evident that at longer cycle lengths, late repolarization was faster in the case without INCX, with APD90 at 0.5 Hz decreasing from 22 to 15 ms, whereas at shorter cycle lengths, repolarization dynamics were less affected by the presence of INCX. The frequency-dependent prolongation in APD70 and APD90 with increasing cycle length was much attenuated in the absence of INCX. Therefore, our model results matched well with the finding that INCX plays a role in shaping the later repolarization phase of the AP.
The electrical restitution of the AP is shown in Fig. 9, where S1 is for the last trace of a pacing train at a cycle length of 300 ms followed by a premature beat (S2), with various S1-S2 intervals ranging between 70 and 1,000 ms. Figure 9B shows the relationship between APD90 and the S1-S2 interval, with APD90 gradually lengthening from 14.8 ms at a 70-ms S1-S2 interval to 16.5 ms at a 1,000-ms S1-S2 interval. This was also in qualitative agreement with experimental observations reported by Knollmann et al. (30) in the intact murine heart. The effects of a premature beat on ICaL and [Ca2+]i is shown in Fig. 9. The APs during the last cycle of the pacing train followed by a premature beat 100 ms later are shown in Fig. 9C, where there was a slight reduction in overshoot during S1 compared with S2 (30.5 vs. 35.8 mV), a slower upstroke (Vmax: 132 vs. 168 V/s), and a faster repolarization with APD70 decreasing from 5.70 ms during S1 to 5.65 ms during S2 and APD90 from 14.9 to 14.4 ms. The changes in the AP shape was mirrored by a decrease in the magnitude of ICaL from 12.0 pA/pF during S1 to 9.8 pA/pF during S2, as shown in Fig. 9D, top. The reduction in Ca2+, in turn, triggered a smaller amount of Ca2+ release from the SR; hence, the smaller [Ca2+]i transient shown in Fig. 9D, bottom.
Ca2+ Handling in Mice Overexpressing Canine Cardiac NCX
I-V relationship of INCX.
[Na+]o, extracellular K+ concentration ([K+]o), and [Ca2+]o were adjusted to 140, 6, and 1 mM, respectively, to match the experimental setup of Terracciano et al. (52). [Na+]i (9.6 mM) and [Ca2+]i (0.135 μM in WT and 0.159 μM in TG) were fixed to the same values as those observed experimentally. The voltage-clamp protocol and simulated INCX in the WT and TG models are shown in Fig. 10A. The simulated I-V relationships are shown in Fig. 10B. At positive potentials, simulated INCX in the TG model was, on average, 2.7 times the magnitude of the WT model. At negative potentials, the differences in INCX were less significant, with an average ratio of 1.5, consistent with experimental observations showing that the statistically significant difference in INCX diminishes toward more negative potentials (52).
Field-stimulated and caffeine-induced [Ca2+]i transients.
Experimentally measured [Ca2+]i transients, as shown in Table 5, showed quite significant variations across different groups resulting from, at least partially, differences in the experimental setup. For example, experiments have been carried out at room temperature (1, 52, 59) or 37°C (53); [Ca2+]i transients were elicited by either depolarizing pulses or field stimulation at frequencies ranging from 0.25 Hz (59) to 0.5 Hz (52) to 1 Hz (53), and perfusion solutions with different [Na+]o, [K+]o, and [Ca2+]o have been used. However, ventricular myocytes from NCX-overexpressing mice have been consistently characterized by a two- to threefold increase in the forward function of NCX (52, 53, 59), an unchanged (1, 52, 53) or a slightly higher (59) [Ca2+]i transient magnitude, a faster rate of decay of the [Ca2+]i transient elicited by either field stimulation or caffeine, and unchanged (1, 59) or enhanced (52, 53) SR Ca2+ content.
To validate our models, the protocol described by Terracciano et al. (52) was reproduced in both the WT and TG models. Results were compared with those reported by Terracciano et al. (52) as well as other studies (1, 53, 59). The extracelluar ion concentrations in both the WT and TG models were also adjusted as previously described. The models were initially paced at 0.5 Hz until the limit cycle was reached. Two seconds after the cessation of stimulation, the tranmembrane voltage was clamped at −75 mV, and the effects of caffeine application were simulated by setting the open probability of RyR to 1 and increasing the rate of SR Ca2+ leak to 0.1 ms−1, thus causing a rapid and complete release of Ca2+ from the SR.
Simulated [Ca2+]i transients during 0.5-Hz field stimulations are shown in Fig. 11A. The magnitude of the [Ca2+]i transient was 0.28 μM in the TG model, 22% higher compared with the WT model (0.23 μM). In general, the [Ca2+]i transient by field stimulation in the TG mouse has been found to be of a similar magnitude or with a slight increasing trend compared with that in WT mice (see Table 3). Simulated peak forward-mode INCX (Fig. 11C) during the [Ca2+]i transient was 0.81 pA/pF in the TG model compared with 0.33 pA/pF in the WT model, corresponding to a 2.5-fold increase in INCX magnitude. The decay of the [Ca2+]i transient was faster in the TG model, with a 30% decrease in RT50 (92 vs. 125 ms) and τ (168 vs. 220 ms). Terracciano et al. (52) reported a 38% decrease in RT50. In another study (53) of the TG mouse, a 33% decrease in RT50 and a 23% decrease in τ have been reported.
Simulated caffeine-induced [Ca2+]i transients are shown in Fig. 11B. The magnitude of the [Ca2+]i transient in the TG model was 0.40 μM compared with 0.30 μM in the WT model, indicating a greater SR Ca2+ load in the TG model. This was in qualitative agreement with the 69% increase in caffeine-induced [Ca2+]i transient magnitude reported by Terracciano et al. (52). In other studies, one study (53) reported a 77% increase, whereas others (1, 59) observed no statistically significant difference in caffeine-induced [Ca2+]i transient magnitude. Simulated peak INCX (Fig. 11D) during caffeine application was 1.44 pA/pF in the TG model compared with 0.44 pA/pF in the WT model, corresponding to a threefold increase in INCX magnitude. A threefold increase in INCX magnitude during caffeine application was reported by Terracciano et al. (52) and two other groups (1, 59), whereas a 2.5-fold increase was reported in another study (53). The rate of decay of the simulated caffeine-induced [Ca2+]i transient was 61% faster in the TG model (τ: 402 vs. 1,022 ms). Although not quantified by Terracciano et al. (52), τ is reportedly 73% and 60% faster in TG mice according to Yao et al. (59) and Adachi-akahane et al. (1), respectively.
Effects of allosteric regulation.
To determine whether the presence of allosteric regulation in the canine cardiac NCX, which is absent in WT mice, plays a significant role in Ca2+ handling in the TG mice, we replaced the current formulation of NCX in the TG model with an alternative formulation (TG*) without allosteric regulation, i.e., VNCXmax was set to 9.75 pA/pF to represent a 2.5-fold increase in maximum NCX activity, and all other parameters values were the same as those in the WT model. In contrast with previous observations, field-stimulated [Ca2+]i transient magnitude (Fig. 11A) was significantly lower in the TG* model compared with the WT model (0.06 vs. 0.23 μM). Caffeine-induced [Ca2+]i transient magnitude (Fig. 11B) was only 23% of the level in WT mice (0.07 vs. 0.30 μM), indicating a dramatic loss of SR Ca2+. Despite the increase in VNCXmax, inward INCX magnitude during field stimulation (Fig. 11C) was lower in the TG* model during systole (0.25 vs. 0.33 pA/pF) and slightly higher during diastole (0.052 vs. 0.043 pA/pF). Consistent with the lower caffeine-induced [Ca2+]i transient magnitude, INCX during caffeine application (Fig. 11D) was also lower in the TG* model (0.30 vs. 0.44 pA/pF), reflecting a significant reduction in SR Ca2+ content. When the rates of decay of the [Ca2+]i transient were compared, the time constant in the TG* model was longer during field stimulation (324 vs. 220 ms) and shorter during caffeine application (628 vs. 1,022 ms).
Species- and temperature-dependent differences in Ca2+ dynamics are well known. However, previous modeling studies have often had the limitations of using heterogeneous experimental data to develop model components (Fig. 1), which precludes the model's potential to reliably interpret experimental observations. The first aim of this study, therefore, was to use consistent experimental data, whenever possible, to parameterize model components and provide a transparent link between parameters and experimental data. All the measurements on [Ca2+]i transients and ICaL were obtained in our own laboratory from LV ventricular myocytes obtained at 35°C. Other transmembrane currents and AP measurements were taken from existing literature, from murine ventricular myocytes at a temperature between 32 and 35°C. Nevertheless, it is important to address the following limitations and discuss the scope in which our model can be used for the interpretation of experimental observations.
First, measurements of the outward K+ currents obtained at physiological temperatures are less abundant than those at room temperature. In some cases, this necessitates the use of data from multiple mouse strains. These different genetic backgrounds are known to have different AP characteristics (9, 55). These differences are likely to be less than differences observed between species but still need to be considered. This model is thus only an approximation of the C57BL/6 murine cardiac myocyte based on currently available data and will need further updating as new data from C57Bl/6 LV ventricular myocytes become available.
Second, we used a model of intracellular Ca2+ buffering proposed by Trafford et al. (54) based on experimental data from ferret ventricular myocytes at room temperature. Whether there are significant differences in buffering properties between mice and ferrets is unclear. To estimate the differences, if any, between ferret and murine Ca2+ buffering, we used simultaneous measurements of caffeine-induced [Ca2+]i transients and INCX in murine ventricular myocytes at room temperature performed by Adachi-Akahane et al. (1) to estimate the murine cardiac myocyte Ca2+-buffering curve. The buffering curve obtained from the Adachi-Akahane et al. (1) data was in close agreement with the Trafford and Eisner buffering curve for [Ca2+]i below 0.5 μM and was slightly lower at higher [Ca2+]i levels. Since in most cases the observed [Ca2+]i transient in our sample had a peak at or below 0.5 μM, the use of ferret buffering parameters seems a reasonable approximation until specific buffering measurements in murine hearts at physiological temperatures become available.
Third, apart from being involved in frequency-dependent accelaration of relaxation, CaMKII is thought to regulate several other Ca2+-transport proteins, including RyRs and the L-type Ca2+ channels. Most studies on CaMKII have been done using acute CaMKII overexpression and endogenous or exogenous activation of CaMKII. For example, Maier et al. (38) found a fourfold increase in resting Ca2+ leak from the SR and a decrease in SR Ca2+ content in TG mice overexpressing CaMKII compared with WT littermates. In another study (19), endogenous CaMKII activation in mice was found to increase Ca2+ spark frequency in both WT and phospholamban knockout mice. Although the exact level of interactions between CaMKII and RyRs under physiological conditions is not clear and was not directly observed in our experiments, we tested the sensitivity of our model to increases in the SR Ca2+ leak rate and the rate of SR Ca2+ release. Doubling the rate of SR Ca2+ release from 4.5 to 9 ms−1 led to no change in the magnitude of the transient at 3 Hz but a 19% decrease in the time to peak and a 10% decrease SR Ca2+ content. This is consistent with our results (described above) showing that changes in the open probability of RyRs only have a transient effect on the size of the [Ca2+]i transient. On the other hand, doubling the rate of SR Ca2+ leak from 3 × 10−5 to 6 × 10−5 ms−1 resulted in a 13% decrease in the magnitude of the transient, a 20% decrease in the SR Ca2+ content, and no change in the time to peak.
The modulation of L-type Ca2+ channels by CaMKII has been reported by several groups. Generally, CaMKII is thought to facilitate ICaL with repeated depolarizations [for a review, see Maier and Bers (37)]. More recently, ICaL in ventricular myocytes from mice with CaMKII overexpression have been reported to have a higher peak and slower decline (38). In rabbit ventricular myocytes, acute CaMKII overexpression led to an increase in peak ICaL and a prolongation of the inactivation kinetics without causing statistically significant changes in the [Ca2+]i transient (31). At 3 Hz, increased activity of CaMKII should be expected to facilitate ICaL; however, both Dibb et al. (13) and Antoons et al. (2) have reported progressive reductions in peak ICaL with increasing frequencies in rat and mouse ventricular myocytes, respectively. Thus, since the physiological role of ICaL facilitation is still debatable (37), it was not included into our model formulation.
Finally, care should be taken when extrapolating Ca2+-handling properties observed in isolated myocytes to in vivo functions. As mentioned in materials and methods, normal contractions at a physiological pacing frequency (10 Hz) could not be sustained in isolated myocytes; thus, we could only obtain [Ca2+]i transients at lower pacing frequencies (up to 3 Hz). This could be due to a number of factors, including cell damage during isolation, loss of structural integrity, and compromised metabolic properties (48, 50). It is therefore possible that other aspects of Ca2+-handling properties are also affected.
Despite the above limitations, simulated Ca2+ dynamics showed good agreement with experimental observations. The description of the RyR in our model was consistent with experimental observations. First, the graded response of Ca2+ release could be reproduced using our model, which showed a bell-shaped voltage relationship with both peak ICaL and the peak Ca2+ transient occurring at ∼0 mV. Second, increasing the open probability of the RyR, by means of caffeine, for example, produced a transient increase in [Ca2+]i transient, which returned gradually to the control level, similar to that observed experimentally by Eisner et al. (15). Net sarcolemmal flux and SR Ca2+ content also followed the same changes as those described by Eisner et al. (15), i.e., a transient net sarcolemmal Ca2+ efflux led to a sustained loss of SR Ca2+ content, which reset the balance point of Ca2+ dynamics after caffeine application.
To achieve temperature consistency across different model components, we also refitted the K+ and Na+ channels to experimental data at physiological temperatures close to 35°C. Experimentally measured AP shape and APD were used to validate the fitting process. There was good agreement in AP characteristics between simulation and experimental data, including resting potential, overshoot, maximum upstroke velocity, and APD25, APD50, and APD90. With increasing cycle length, steady-state APD increases were in agreement with observations made by other groups (3, 30, 51). This was shown to be strongly affected by the presence of INCX in our model, which was also seen experimentally (28). Therefore, our model is capable of representing the interactions and feedback loops between Ca2+ dynamics and sarcolemmal electrophysiology at physiological pacing frequencies and temperatures under steady-state conditions.
The second aim of this study was to use our model to investigate the effects of heterozygous overexpression of canine NCX in TG mice. It is worth noting that the more recent findings by Weber et al. (56) on allosteric regulation are not fully consistent with those reported previously by Maxwell et al. (39), who observed allosteric [Ca2+]i regulation of NCX in both control and TG mice. This indicates that the lack of [Ca2+]i regulation observed by Weber et al. (56) in WT preparations may be due to species- and/or preparation-dependent differences. However, as the study by Weber et al. (56) was conducted under more physiological conditions (in intact myocytes and allowing dynamic [Ca2+]i regulation) compared with previous measurements obtained in membrane patches, we parameterized the present models based on their measurements.
Our simulation results from the TG model were consistent with the experimentally observed 2.7-fold increase in NCX activity (1, 52, 53, 59) assessed from the I-V relationship and the rate of decay of the caffeine-induced [Ca2+]i transient. At the pacing frequency of 0.5 Hz, the simulated [Ca2+]i transient magnitude in the TG model was 22% greater compared with the WT model with a 30% faster rate of [Ca2+]i decay and a 33% increase in SR Ca2+ load. These results were in agreement with experimental studies (1, 52, 53, 59).
No difference in the time course of ICaL was observed between the WT and TG models stimulated at 0.5 Hz (not shown here), so the total amounts of Ca2+ entry through the sarcolemma are the same between WT and TG models. The balance of Ca2+ entry and extrusion over each heart beat, therefore, require that Ca2+ extrusion through NCX should also be the same in the two models. Indeed, our calculations of the integral of JNCX showed the same amounts of Ca2+ extrusion. This occurred despite the over twofold increase in the magnitude of JNCX in the TG model, which was compensated by a lower diastolic JNCX. The reason behind the decrease in JNCX during diastole was not a lower diastolic [Ca2+]i, as diastolic [Ca2+]i was, in fact, slightly higher in the TG model at 0.5 Hz. However, this may be explained by the presence of allosteric regulation by [Ca2+]i in TG mice, which has not been found in WT mice (56). According to the NCX formulation (see Supplemental Material), with an estimated KmAllo of 0.15 μM (56) and a diastolic [Ca2+]i of 0.05 μM in our models, the maximum exchange rate of NCX could be as low as a 10th of VNCXmax during diastole in the TG model. However, during systole, when [Ca2+]i is high, the maximum exchange rate approaches VNCXmax, so peak INCX is considerably higher in the TG model.
The less competitive diastolic JNCX and higher diastolic [Ca2+]i promote Ca2+ uptake into the SR through SERCA. This means that more Ca2+ is available for CICR and, hence, there is a higher [Ca2+]i transient magnitude in the TG model. The greater [Ca2+]i transient magnitude, in turn, facilitates SR Ca2+ uptake, the consequences of which are twofold. First, higher JSERCA, along with the enhanced systolic JNCX, results in a faster decay of the [Ca2+]i transient. Second, it increases the relative contribution of SERCA and, hence, decreases that of NCX in the TG model. Although a difference in diastolic [Ca2+]i has not been detected experimentally (1, 52, 59), Terracciano et al. (52) speculated that a slightly higher diastolic [Ca2+]i would promote SR uptake and SR Ca2+ loading in TG mice, although the difference might be smaller than can be measured.
In contrast with the increased [Ca2+]i transient magnitude, the increased SR Ca2+ load, and the faster decay kinetics of [Ca2+]i observed in the TG model, results from the TG* model show different characteristics of Ca2+ handling. When allosteric regulation was removed, the enhanced NCX activity alone caused significantly decreased [Ca2+]i transient magnitude and SR Ca2+ content compared with the WT model. This occurs as JNCX is enhanced during both diastole and systole in the absence of allosteric regulation, therefore shifting the balance away from SR Ca2+ uptake and toward sarcolemmal Ca2+ extrusion. This, in turn, leads to a net loss of Ca2+ and lowered systolic [Ca2+]i to a point where sarcolemmal Ca2+ influx and efflux are again in balance. The decrease in [Ca2+]i transient magnitude also decreases JSERCA such that the rate of decay of the [Ca2+]i transient is slightly slower in the TG* model during 0.5-Hz field stimulation. The above simulation results suggest the presence of allosteric regulation in the canine cardiac NCX could be a key factor underlying the observed increases in [Ca2+]i transient magnitude and SR Ca2+ content in TG mice.
The main contribution of the present study is the development of a species-specific and temperature-consistent model framework of Ca2+ dynamics in murine ventricular myocytes. To achieve this goal, available consistent experimental data sources from the mouse were applied within a parameterization methodology developed to establish a transparent link between data and parameter values. The model was validated against a range of experimental measurements, including the frequency response of Ca2+ dynamics, AP shape, the dependence of APD on cycle length, and electrical restitution. The present framework allows the characterization of key Ca2+-handling parameters using a set of data from standard experimental procedures. The model can also be used in many applications, such as studies involving TG mouse models, to allow quantitative analysis of the function. As a case study, TG mice overexpressing the canine NCX isoform were investigated using our model. Simulation results indicated that observations from genetic modifications can be species dependent, and in the case where species variation in the NCX isoform is eliminated, acute overexpression of NCX may lead to a net loss of Ca2+ from the cell, which lowers both SR Ca2+ load and systolic [Ca2+]i transient magnitude. This result further supports the need to use a species-specific model that takes into account species-dependent variations to help interpret experimental observations and elucidate the underlying mechanisms that may be difficult to observe directly.
No conflicts of interest, financial or otherwise, are declared by the author(s).
L. Li thanks the Clarendon Fund for providing funding for this research and Dr. Jonathan Cooper for help with the sensitivity analysis. S. A. Niederer acknowledges the support of the Engineering and Physical Sciences Research Council. N. P. Sweitach and B. Casadei acknowledge the support of the Medical Research Council.
↵1 Supplemental Material for this article is available online at the American Journal of Physiology-Heart and Circulatory Physiology website.
- Copyright © 2010 the American Physiological Society