Integrative Systems Models of Cardiac Excitation–Contraction Coupling

Common Pool Models of CICR

The landmark DiFrancesco–Noble model of the cardiac Purkinje fiber²⁸ was developed from the preceding Beeler–Reuter model²⁹ (first computational model of the cardiac ventricular myocyte) and that of McAllister et al.³⁰ This model was ground breaking for many reasons and was the first model to introduce descriptions of the transport of Ca²⁺ into the NSR by SERCA, extrusion of Ca²⁺ from the cytosol by NCX, diffusion of Ca²⁺ from NSR to JSR, and CICR from the JSR. This model established the conceptual framework on which all subsequent models of the myocyte have been built. Although this model did reproduce graded SR Ca²⁺ release, it did so with ECC gain that was nonphysiologically low.

The construction of mechanistic mathematical models of CICR which attempt to capture both high gain and graded SR Ca²⁺ release first lead to a series of “common pool” models, which, in fact, fail to capture these fundamental properties. As defined by Stern,²⁷ common pool models are those in which LCC trigger Ca²⁺ enters the same Ca²⁺ pool into which JSR Ca²⁺ is released, and which controls JSR Ca²⁺ release. The rapid increase of Ca²⁺ in this pool leads to nonphysiological regenerative, all-or-none rather than graded Ca²⁺ release, a situation where whole-cell SR Ca²⁺ release escapes control by I_CaL once it is initiated.²⁶ An additional assumption inherent in this type of model is that the gating dynamics of the entire population of RyRs are controlled by properties of whole-cell I_CaL rather than by properties of local LCC channel gating. Stern elegantly demonstrated via the application of linear stability theory that common pool models cannot achieve both high gain and graded Ca²⁺ release.²⁷

Regardless of this limitation of common pool models, many advances in cardiac myocyte modeling over the past decade have stemmed from improved descriptions of CICR in such models. The Jafri et al³¹ model of the guinea pig ventricular myocyte was the first to incorporate a Markov model of the LCC (Figure 1A) into an integrated myocyte model with mechanistic ECC. The LCC model was based on the concept of modal gating introduced by Imredy and Yue³² based on single-LCC recordings in which Ca²⁺ induced a shift in LCC-gating dynamics into a low open probability mode. The model contains subsets of states which reproduce channel behavior before (mode-normal) and following (mode-Ca) the occurrence of Ca²⁺ dependent inactivation (CDI). Openings are rare in mode-Ca, effectively inactivating the channel. Successive LCC models have incorporated the concept of modal gating in Markov models for describing LCC inactivation^33,–,35 as well as other modes such as phosphorylation.³⁶ The Jafri et al³¹ model was also first to implement a “restricted subspace,” a single compartment representing the total volume of all dyads into which all Ca²⁺ flux through RyRs and LCCs is directed. This model made a number of predictions regarding mechanisms of CICR and Ca²⁺ cycling (reviewed elsewhere³⁷). A major prediction of this model was that experimentally measured interval-force (transient) and frequency-force (steady-state) relationships could be explained by the interplay between RyR inactivation and SR Ca²⁺ load dynamics. This Ca²⁺ subsystem was modified for canine in the model of Winslow et al³⁸ from which a model of heart failure was developed based on data obtained from myocytes isolated from dogs with tachycardia pacing–induced heart failure.³⁹ The model predicted that in heart failure, downregulation of the SERCA pump and upregulation of NCX decreased JSR Ca²⁺ load, Ca²⁺ release via the RyRs, and the strength of I_CaL CDI. The resulting increase in I_CaL acted to prolong AP duration (APD), a characteristic phenotype of heart failure. The model predictions were subsequently validated by experimental studies,^40,41 thereby demonstrating that negative feedback on LCCs by Ca²⁺ release from RyRs can exert significant control of AP properties. Furthermore, the heart failure model could produce early afterdepolarizations (EADs), a cellular arrhythmia believed to be a possible trigger for arrhythmia at the whole heart level. On introduction into a model of the ventricles, the EADs in the heart failure model quickly led to degeneration of electric propagation into a cardiac arrhythmia, whereas the control myocyte model yielded normal coordinated ventricular activation during pacing.⁴²

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

A, Markov state model for gating of the LCC as proposed by Jafri et al.³¹ CDI occurs as a result of downward transitions from mode-normal to mode-Ca. VDI is described by an independent Hodgkin–Huxley type gate (not shown). B, Structure of a region of the dyad containing 1 LCC opposed to a single RyR, and 1 CaM molecule tethered to the LCC.⁷³ Reprinted with permission from Elsevier.

Despite the fact that these common pool models were not able to reproduce graded SR Ca²⁺ release, they were able to reproduce and predict many important behaviors of the cardiac ventricular myocyte, raising the question of whether capturing the mechanisms of graded release is of primary importance. One way to circumvent this issue is to model SR Ca²⁺ release flux as a phenomenological function of LCC influx and/or membrane potential. Doing so removes the positive feedback effect inherent to common pool models. Models incorporating a phenomenological mechanism of CICR have been successful in describing many myocyte behaviors (including the widely used models of Rudy and colleagues^43,44). However, as with any phenomenological model, care must be taken when interpreting results of these models that pertain to Ca²⁺ release. The critical relationship between SR Ca²⁺ release and APD,^38,40,41 and the recognition that the feedback of SR Ca²⁺ release on LCCs occurs locally in discrete Ca²⁺ spark events²¹ suggested that capturing the mechanisms of local interaction between LCCs and RyRs would be essential. Furthermore, experimental data were emerging that supported the idea that LCC inactivation was dominated by CDI rather than voltage-dependent inactivation (VDI).^45,–,47 Strong CDI in combination with weak/incomplete VDI causes instability in common pool models³⁴ (which were therefore typically formulated with strong VDI and weak CDI) because all-or-none SR Ca²⁺ release induces fast whole-cell CDI which eliminates the late phase of I_CaL and destabilizes the AP plateau (see figure 10C in the article by Greenstein and Winslow³⁴). This made it clear that to fully capture the core mechanisms of ECC, a model in which LCCs locally control RyR function was needed.

Local Control of SR Ca²⁺ Release

In 1992, Stern²⁷ formulated idealized models of “local control” of CICR, and in doing so advanced the local control theory of ECC which underlies how the microscopic properties of CICR translate into macroscopic physiological phenomena observed in experiments. The approach was to recognize that the nonphysiologic positive-feedback loop inherent in common pool models of CICR must be broken by separating Ca²⁺ released from the SR from the trigger Ca²⁺ influx. This was accomplished by noting that LCCs and RyRs reside in close apposition on opposite sides of the dyadic cleft, which would effectively grant Ca²⁺ ions permeating an LCC privileged access to the Ca²⁺ binding sites on nearby RyRs. Importantly, the assumed structural separation of individual clefts⁸ means that released Ca²⁺ would not generally act as a trigger to RyRs in distant dyads (although recent studies suggest neighboring RyR clusters may have a functional role in ECC^11,48). Stern's “calcium-synapse” model, in which each dyad was assumed to contain a single LCC and a single SR release channel, captured these ideas in an elegantly simple way, demonstrating that high gain and graded SR release arise from the local correlation in stochastic LCC and RyR gating because of their communication via a common local Ca²⁺ signal. In this model, graded release is achieved by statistical recruitment of elementary SR Ca²⁺ release events. This fundamental idea of local control of SR release spurred a number of modeling studies, some of which focused on detailed dynamics of Ca²⁺ distribution within the dyadic junction, some of which focused on mechanisms of Ca²⁺ release termination, and a small handful of which developed fully integrative models that relate macroscopic cellular properties of ECC to the function of independent, locally controlled Ca²⁺ release units.

Because of the small size of the dyad, it became clear that understanding the spatial and temporal profile of [Ca²⁺]_d would be an important step in elucidating the mechanistic details of local CICR. Langer and Peskoff⁴⁹ developed one of the first detailed models to numerically simulate the time course of dyadic Ca²⁺. The model predicted that a single typical LCC opening led to ≈1 mmol/L [Ca²⁺]_d near the mouth of the LCC. The presence of sarcolemmal buffering sites slowed the decay of the local Ca²⁺ transient (such that it remained highly elevated long after LCC closure) as compared to the dynamics in the absence of buffers, a finding that was inconsistent with the idea that LCC gating imposes tight control on RyR activity. Soeller and Cannell⁵⁰ formulated a dyad model that introduced the electrostatic effects on ion movement that predicted that [Ca²⁺]_d actually changes more rapidly in response to LCC gating. In this model, [Ca²⁺]_d was approximately proportional to unitary LCC currents over a wide range of current amplitudes. The model yielded predictions of [Ca²⁺]_d profiles during JSR release, which peaked as high as 300 μmol/L near the center of the dyad. The decay in dyadic Ca²⁺ occurred in ≈2 ms following termination of SR release (but this was sensitive to the dyad boundary condition). The high speed of the response of dyadic Ca²⁺ to LCC and RyR gating in this model helped to explain how tight local control of SR release by LCCs can be achieved.

With the growing acceptance of local control theory of CICR and that Ca²⁺ sparks are the elementary events underlying SR Ca²⁺ release (for a review see⁵¹), mathematical models of Ca²⁺ sparks were developed. The summation of Ca²⁺ sparks could explain the characteristic properties of CICR, ECC gain, and whole-cell Ca²⁺ transients; therefore, the analysis of models of single Ca²⁺ sparks and ensembles of Ca²⁺ sparks provided a wealth of information. Local control of CICR requires that SR release be stable, which suggests that Ca²⁺ sparks, which are locally regenerative, must be terminated rapidly by a robust mechanism. Although the mechanisms of Ca²⁺ spark termination still remain unclear, a number of models have helped elucidate how possible mechanisms may contribute to termination of SR Ca²⁺ release.

Stern²⁷ originally described 3 mechanisms that can terminate SR release in a cluster of RyRs: (1) RyR inactivation is a time-dependent process which may be initiated by Ca²⁺ binding to an inactivation site on the RyR or may occur “fatefully” following opening of the channel. The observed inactivation process has been described as classically inactivating (leading to absolute refractoriness of CICR)⁵² or taking on the complex properties of adaptation.⁵¹ (2) Depletion of locally releasable JSR Ca²⁺ will occur as the result of RyR activity. The emptying of the cisternal Ca²⁺ store reduces the Ca²⁺ gradient across the RyRs and hence the Ca²⁺ flux. This will diminish the local feedback that underlies regenerative release in a cluster of RyRs, allowing them to close. (3) “Stochastic attrition” occurs as a result of the nonzero probability that all RyRs in a cluster are simultaneously closed. Because Ca²⁺ diffuses rapidly out of the dyad, the simultaneous closure of all RyRs would lead to a rapid reduction in [Ca²⁺]_d below the level which sustains RyR activity. Active RyR clusters would therefore gradually cease gating spontaneously. This mechanism can only be effective on a reasonable timescale if very few RyRs are activated during a Ca²⁺ spark, and it does not provide for a mechanism of refractoriness. The stereotypical nature of Ca²⁺ spark properties suggested that spark generation and termination are likely controlled by regulatory feedback mechanisms and cannot be described by a simple stochastic mechanism alone.²³ An additional mechanism proposed recently is RyR regulation by luminal SR Ca²⁺. RyR sensitivity to [Ca²⁺]_d can be altered by the level of JSR luminal [Ca²⁺] via the interaction of the Ca²⁺ buffer calsequestrin with RyR accessory proteins triadin and junction within the JSR.⁵³ The decrease in RyR open probability with JSR depletion likely contributes to termination of Ca²⁺ release.

In 1999, Stern et al³³ developed a stochastic model of cardiac ECC. The model consisted of a cardiac dyad in which the Ca²⁺ concentration profile was simulated by solving equations governing diffusion and Ca²⁺ binding reactions in the cleft. LCCs and RyRs were represented by Markov state models, which were simulated as part of a stochastic system in which the state of each individual channel was tracked. The behavior of 6 different RyR gating schemes was analyzed in the context of local control of ECC. This study revealed that local control, in a realistic setting of dyad geometry, number, and location of channels, could explain graded release, when the RyR was represented by a simple phenomenological four-state gating scheme which included inactivation (Ca²⁺-dependent or fateful). This success has led to continued development of this RyR model in more recent work.^54,55 However, RyR channel gating schemes that were based on single-channel RyR measurements in lipid bilayers^56,57 failed to give stable ECC (but stability could be restored by introducing allosteric RyR interactions³³). The work of Stern et al³³ suggested that there must be a sufficiently strong process that substantially reduces RyR open probability shortly following activation such that stochastic attrition can occur on the timescale of physiological ECC.

Whether or not RyR inactivation or adaptation plays an important role in termination of SR Ca²⁺ release remains controversial (see below). A number of RyR models^57,58 were developed to better understand the observation of RyR adaptation,⁵⁹ whereby an RyR can reopen in response successive increases in Ca²⁺ concentration, indicating they have not entered a conventional absorbing CDI state. Rice et al⁶⁰ developed a stochastic model of Ca²⁺ release in the dyadic junction which demonstrated that RyR adaptation modulated Ca²⁺ release duration and amplitude but was not required for termination of Ca²⁺ release (ie, JSR depletion played an important role). This model achieved high gain and graded release without tracking the spatial gradients of [Ca²⁺]_d. In contrast, the RyR model of Shannon et al⁵⁴ is a modification of the phenomenological four-state model of Stern et al,³³ which introduced an important modulatory role for SR luminal Ca²⁺ in RyR gating but which still depends on an inactivation/adaptation mechanism for release termination. Variations of this model have been incorporated into a variety of recent models which have successfully described SR load and frequency dependent properties of ECC and APs in rabbit⁵⁴ and human,⁶¹ the role of RyR phosphorylation on Ca²⁺ dynamics and the AP,⁵⁵ and the role of spatial coupling between dyadic release sites in Ca²⁺ spark–induced sparks and Ca²⁺ alternans.⁴⁸

In 2002, Sobie et al⁶² developed a Ca²⁺ spark model to further investigate mechanisms underlying termination of SR Ca²⁺ release. In this model, each release site is composed of a large number of two-state RyRs which incorporate 2 additional experimental observations: (1) elevated SR luminal Ca²⁺ increases the Ca²⁺-sensitivity of RyR activation⁶³; and (2) RyR coupling cooperatively influences their gating properties within a cluster.⁶⁴ These mechanisms of RyR gating yield a “sticky cluster” model which reproduces fundamental features of Ca²⁺ sparks and robust termination of release that arises as a result of the combined effects of these 2 mechanisms. In this model, the assumption of coupled RyR gating is essential for reliable SR release termination. However, it remains unclear how RyRs would influence the gating properties of all other RyRs in a cluster (as opposed to only adjacent RyRs³³) as implemented in this model (see elsewhere^65,66 for alternative approaches to modeling the functional effects of RyR coupling). In addition, the >90% SR depletion required for release termination in the Sobie et al⁶² model may not be realistic as experiments have demonstrated that 40% to 75% of releasable Ca²⁺ remains in the JSR following whole-cell triggered SR release.^67,68 Despite these issues, the majority of recent experiments^68,–,70 have strengthened the idea that termination of SR release is primarily controlled by luminal SR Ca²⁺ rather than an RyR inactivation mechanism.

Although many controversies as to the details of CICR remain to be resolved, the models of dyadic Ca²⁺ dynamics, Ca²⁺ sparks and CICR have helped affirm or disprove a variety of mechanistic hypotheses regarding cardiac ECC. Most of these ECC models, however, have been implemented at the subcellular scale, isolated from the physiological environment and influence of other cellular processes. The first mathematical model to incorporate the mechanisms of local control of SR Ca²⁺ release at the level of the cardiac dyad into an integrative model of the cardiac action potential was published in 2002 by Greenstein and Winslow³⁴ (Greenstein–Winslow model). This model is described in detail in the following section.

Multiscale Modeling of CICR

Quantitative understanding of CICR requires that Ca²⁺ signaling be understood across a wide range of spatiotemporal scales. The length scale over which CICR occurs within the dyad is on the order of nanometers, and relevant time scales are of the order of microseconds. Integration to the level of the cell requires formulation of models that describe behavior at spatiotemporal scales of several hundred microns and tens to hundreds of seconds. Recently, we have developed several new approaches for multiscale modeling of CICR spanning these spatiotemporal scales.^34,71,–,73 In this section, we present these models and the mathematical techniques for moving across these spatiotemporal scales.

The Nanoscale Model of CICR

We have recently developed a computational model of Ca²⁺ signaling in the dyad that simulates motion of Ca²⁺ ions, binding of ions to RyR binding sites, and Ca²⁺ release from RyR. This model was developed to understand factors shaping CICR at the nano-scale level. We addressed the following 3 questions. (1) What is the number of Ca²⁺ ions that are present in the dyad during CICR? (2) How does the physical arrangement of large proteins within the dyad influence CICR? (3) How does “signaling noise” caused by the potentially small number of Ca²⁺ ions within the dyad affect the nature of CICR? We addressed these questions by developing a molecularly and structurally detailed computational model of the cardiac dyad describing motions of Ca²⁺ ions, as influenced by thermal energy and electrostatic potentials.

The small size of the dyadic cleft (≈4×10⁵ nm³ with ≈10 to 40 RyRs),^10,11 in combination with the relatively large proteins that reside in the cleft, results in a space in which Ca²⁺ diffusion is highly restricted. Figure 1B shows a representation of an LCC, RyR, and a single calmodulin (CaM) molecule in the dyad based directly on structural data.^74,–,76 Models predict that during an AP, peak Ca²⁺ concentration in the dyad ranges from 100 to 1000 μmol/L,^49,50 which corresponds to as few as 20 to 200 free Ca²⁺ ions.⁷³ Therefore, application of laws of mass-action to analysis of signaling events mediated by this relatively small number of ions may not be justified.⁷⁷ Rather, the stochastic motion of Ca²⁺ ions in the dyad may impart a degree of signaling noise between LCCs and RyRs, thereby influencing properties of CICR.⁷³ In addition, the physical location and dimensions of dyad proteins may have considerable influence on motion of Ca²⁺ ions and hence, CICR.

Development of the nanoscale model is presented in detail elsewhere.⁷³ Briefly, a simplified dyad was modeled using a 200×200×15 nm lattice containing 20 RyRs arranged in a 4×5 matrix. Four LCCs were positioned randomly across this lattice to achieve a 5:1 RyR:LCC ratio. The geometry of each LCC and RyR protein was modeled as shown in Figure 1B, with each LCC tethering with one CaM. LCC gating was simulated stochastically using the model of Figure 1A. The gating kinetics of the RyR were described using the four-state Markov model of Stern.^33,78 Entry of Ca²⁺ ions into the dyad via open LCCs and RyRs was modeled as a Poisson process with rate governed by channel permeabilities measured experimentally.⁷⁹ CDI of an LCC required that the 2 carboxyl-terminal Ca²⁺ binding sites of the associated CaM were occupied. Opening of an RyR required that Ca²⁺ was bound to all 4 activation sites. RyRs were assumed to inactivate when Ca²⁺ was bound to at least 1 of 4 distinct inactivation sites. Sarcolemmal anionic sites acted as Ca²⁺ buffers.^49,50 The joint positions of Ca²⁺ ions present in the dyad were modeled as a Brownian motion in a potential field, describing both interactions of Ca²⁺ ions with other Ca²⁺ ions as well as electrostatic potentials based on the Fokker–Planck equation (FPE).⁸⁰ We generated sample paths of Ca²⁺ ion movement in the dyad using a finite difference representation of the FPE that can be interpreted as a spatially discrete Markov process as described by Wang et al.⁸¹

A 10-ms simulation of a model dyad containing only 1 LCC located at its center at a membrane potential of 0 mV generated an average of ≈0.5 Ca²⁺ ions within a 50 nm radius of the LCC (see figure 8 in the article by Tanskanen et al⁷³), with a minimum of 0 and a maximum of 7 Ca²⁺ ions, corresponding to 98.7 μmol/L [Ca²⁺]_d. A similar protocol using a single RyR in the SR membrane (in which the RyR is initially open) generated an average of ≈22.8 Ca²⁺ ions. The number of Ca²⁺ ions peaked at 49 during Ca²⁺ release, which corresponds to a concentration of ≈0.7 mmol/L Ca²⁺, similar to the concentration of Ca²⁺ found in the JSR.⁸² With a coefficient of variation of 64%, RyR gating produces a high degree of variability in the number of dyadic Ca²⁺ ions, showing that the LCC-RyR signaling involved in CICR at the dyad level is a noisy process mediated by tens of Ca²⁺ ions.

Figure 2A and 2B demonstrates fundamental features of a representative Ca²⁺ release event in a dyad in response to a voltage clamp to 0 mV at time 0 ms. Figure 2A shows the number of Ca²⁺ ions in the dyad, and Figure 2B shows the number of LCCs (gray line) and RyRs (black line) that are open as a function of time during the release event. Two LCCs open at ≈3 ms following the voltage clamp, and Ca²⁺ influx into the dyad triggers the opening of 3 RyRs 1 to 2 ms later, which inactivate after ≈7 ms, ≈9 ms, and ≈13 ms, respectively. The temporal shape of the Ca²⁺ signal in the dyad (Figure 2A) is closely correlated to the number of open RyRs (Figure 2B). Figure 2C shows the average dyadic Ca²⁺ signal, calculated by averaging the number of Ca²⁺ ions at each point in time across 400 independent dyads. The duration of the average local Ca²⁺ transient shown in Figure 2C is ≈20 ms (at half-maximal amplitude), similar to that measured for Ca²⁺ spikes in myocytes.^21,24,52

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

A CICR event in a single dyad evoked by a 0-mV voltage clamp step at time 0 ms (adapted from Tanskanen et al⁷³). A, The number of free Ca²⁺ ions in the dyad volume as a function of time. B, The number of open RyRs (black solid line) and the number of open LCCs (gray solid line) in the dyad during the release event depicted in A. C, The average number of free Ca²⁺ ions in the dyad as calculated from 400 independent dyads. D, ECC gain as a function of membrane potential for the baseline model, which includes space-filling geometric models of protein structure in the dyad (solid line), for the model with protein structures excluded (dashed line), and for a modified model with dyad height reduced from 15 to 13 nm and protein structures excluded (dotted line).

We demonstrated that ECC gain is a robust and reproducible measurement of CICR because the variability in gain that arises from the stochastic gating of LCCs and RyRs is small for a whole-cell population of dyads at potentials near and above 0 mV (see figure 11 in the article by Tanskanen et al⁷³). The functional influence of protein structures on ECC is shown in Figure 2D. Peak ECC gain obtained for the baseline model (solid line) is compared to that for the model in which LCC, RyR, and CaM molecule structures were omitted (dashed line). Ca²⁺ binding site locations for the no-structure simulations remain at the same positions in space as when the protein structures are in place. Over a wide range of potentials, the peak ECC gain is decreased on removal of the protein structures from the dyad. The protein structures occupy ≈15% of the dyad volume. To show that the difference in ECC gain in the presence versus absence of protein structures is not simply attributable to the difference in Ca²⁺-accessible dyad volume, peak ECC gain was obtained with protein structures absent and dyad volume reduced by ≈15% (dotted line). ECC gain values for this model were clearly not increased to the level obtained in the presence of protein structures. These data indicate that the physical shape and configuration of dyad proteins influences Ca²⁺ diffusion during CICR in such a way as to enhance ECC gain. In addition, the probability of triggering a release event was 0.152 with protein structures excluded and 0.105 with protein structures intact, a 44% increase. These data suggest that when LCCs are directly apposed to RyRs, the diffusion obstacle imposed by the structure of RyRs leads to an increase in the probability that trigger Ca²⁺ ions bind to activation sites on the RyRs, and hence increases ECC gain.

These data demonstrate that at the level of single dyads, CICR is a highly noisy process and that models based on laws of mass action may not accurately capture dyadic Ca²⁺ dynamics. Further, in ensembles of dyads, the fundamental properties of local control of CICR arise de novo as a consequence of the underlying physical structure and channel gating properties, including voltage-dependence of ECC gain. Finally, the physical shape and location of RyRs relative to LCCs is likely to have a major influence on properties of CICR.

The Integrated Local-Control Model of the Cardiac Myocyte

The Greenstein–Winslow³⁴ model of the cardiac ventricular myocyte was the first to fully integrate mechanisms of local control of CICR into a model of the cardiac AP. This was accomplished by updating and extending the common-pool model of Winslow et al³⁸ to include a population of dyadic Ca²⁺ release units (CaRUs). In essence, this integrated local control model is generated from the nanoscale model by simplifying the dyad description to that of well-mixed Ca²⁺ compartment, and incorporating a large population of such dyads into a whole-cell model. Detailed properties of the nanoscale model (eg, effect of protein structure on ECC gain) are retained in this integrated model by adjusting parameters that influence the effectively equivalent property at this scale (eg, effective Ca²⁺ sensitivity of RyRs). In this model, local interactions of individual LCCs with nearby RyRs are simulated stochastically, with these local simulations embedded within the numeric integration of the ordinary differential equations (ODEs) describing ionic and membrane pump/exchanger currents, SR Ca²⁺ cycling, and time-varying cytosolic ion concentrations.

The Greenstein–Winslow model was formulated to capture fundamental properties such as graded release, although, at the same time, it needed to be computationally tractable. A full mathematical description of the stochastic state models, dynamical equations, parameters, and initial conditions defining the myocyte model are given in Appendix I of a previously published article.³⁴ Figure 3A shows a schematic of the CaRU model which is intended to mimic the properties of Ca²⁺ sparks in the t-tubule/SR (T-SR) junction. Figure 3B shows a cross-section of the model T-SR cleft, which is divided into 4 individual dyadic subspace compartments arranged on a 2×2 grid. Each subspace compartment contains a single LCC and 5 RyRs in its JSR and sarcolemmal membranes, respectively.⁸³ The division of the CaRU into 4 subunits allows for the possibility that an LCC may trigger Ca²⁺ release in adjacent subspaces (ie, RyR recruitment). On activation of RyRs, subspace [Ca²⁺] will become elevated. This Ca²⁺ freely diffuses to either adjacent subspace compartments (J_iss) or into the cytosol (J_xfer). The local JSR compartment is refilled via passive diffusion of Ca²⁺ from the NSR compartment (J_tr). The number of active LCCs was chosen such that the amplitude of the ensemble current summed over all LCCs corresponds to whole-cell measurements in canine myocytes,⁸⁴ which corresponded to 50 000 LCCs (12 500 CaRUs) and agrees with experimental estimates of LCC density.⁷⁹ LCC and RyR gating were based on previous models,^31,85 and the channels (ClCh) that carry the Ca²⁺-dependent transient outward chloride (Cl^-) current (I_to2) are included within the CaRU.⁸⁶ A detailed description of the local control simulation algorithm is given in Appendix II of a previously published article.³⁴ Stochastic fluctuations in model outputs are the natural result of the ensemble behavior of ion channel and CaRUs, and introduce a degree of variability to simulation output and can be physiologically important to understanding some physiological phenomena such as EADs.⁸⁷

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Schematic representation of the CaRU.³⁴ A, Trigger Ca²⁺ influx through the LCCs enters into the T-SR cleft (dyadic space), RyRs and ClChs open, local Ca²⁺ passively diffuses into the cytosol, and JSR Ca²⁺ is refilled via passive diffusion from the NSR. B, The T-SR cleft (shown in cross-section) is composed of 4 dyadic subspace volumes, arranged on a 2×2 grid, each containing 1 LCC, 1 ClCh, and 5 RyRs. Ca²⁺ in any subspace may diffuse to a neighboring subspace (J_iss) or to the cytosol (J_xfer). Reprinted with permission from Elsevier.

Figure 4A and 4B demonstrate the most elementary model release event during CICR, as triggered by a single LCC at 0 mV. Ca²⁺ flux through an LCC (gray line) and the net SR Ca²⁺ release flux through the 5 adjacent RyRs (black line) are shown in Figure 4A. Local JSR release flux is triggered by the first LCC opening (at ≈5 ms) and lasts ≈20 ms, far longer than the LCC open duration (<1 ms). The amplitude of the release flux varies with the number of open RyRs and the local Ca²⁺ gradient across the JSR membrane. Figure 4B shows the corresponding subspace [Ca²⁺], which reaches a peak value of ≈40 μmol/L. Total Ca²⁺ influx through the set of 4 LCCs (gray line) and the net SR Ca²⁺ release flux through the set of 20 RyRs (black line) that make up a single T-SR junction are shown in Figure 4C. LCC Ca²⁺ influx rises to a level consistent with 2 open channels within a short time following the initiation of the pulse. The spatial average of Ca²⁺ concentration in the 4 subspace compartments of the CaRU is intended to represent a Ca²⁺ spark-like event (Figure 4D), and has duration of 25 ms (at half-maximal amplitude) similar to that measured in myocytes (20 to 50 ms).^21,24,52

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Sample results for a single CaRU in response to a 200-ms voltage clamp to 0 mV.³⁴ A, Ca²⁺ flux through a single LCC (gray line) and through the set of 5 RyRs (black line) within a single dyadic subspace compartment. Arrows 1 and 2 highlight RyR number and Ca²⁺ gradient–driven changes in SR Ca²⁺ release flux, respectively. B, Subspace [Ca²⁺] associated with the events of A. C, Ca²⁺ flux through the set of 4 LCCs (gray line) and the set of 20 RyRs (black line) within a single CaRU. D, Mean subspace [Ca²⁺] in the 4 subspace compartments associated with the events in the CaRU described in C. Reprinted with permission from Elsevier.

Previous experimental studies^25,88 and mathematical models^27,33 have shown significant differences between the voltage dependence of I_CaL and RyR Ca²⁺ release flux even though SR Ca²⁺ release is exclusively controlled by LCC Ca²⁺ entry. These differences underlie the phenomenon of “variable” ECC gain. Figure 5A shows the voltage dependence of LCC (filled circles) and RyR (open circles) Ca²⁺ flux. In Figure 5B, the peak fluxes of Figure 5A have been normalized based on their respective maxima. Maximal LCC Ca²⁺ influx occurs at 10 mV, whereas maximal RyR Ca²⁺ release flux occurs at 0 mV. ECC gain is plotted as a function of voltage in Figure 5C (triangles), and is monotonically decreasing with increasing membrane potential in agreement with experimental measurements.^24,25,88 In the negative voltage range LCC open probability is submaximal, leading to sparse LCC openings (ie, submaximal recruitment of CaRUs), but large amplitude unitary currents efficiently trigger adjacent RyRs, resulting in high values of ECC gain. More CaRUs are recruited with increasing membrane potential as a result of increasing LCC open probability, however ECC gain decreases with increasing membrane potential because the higher LCC open probability leads to redundancy in LCC openings as a result of which a fraction of LCCs are attempting to trigger refractory RyRs and at even higher potentials the decrease in LCC unitary current amplitude reduces the fidelity of local coupling between LCCs and RyRs.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Ca²⁺ fluxes and ECC gain for the Greenstein–Winslow model (A through C)³⁴ and coupled LCC-RyR model (D through F).⁷¹ A, Mean peak LCC (filled circles) and RyR (open circles) Ca²⁺ flux amplitudes as a function of membrane voltage. B, Normalized peak Ca²⁺ fluxes (data of A). C, ECC gain under control conditions. D, Peak LCC (filled triangles) and RyR (open triangles) Ca²⁺ flux amplitudes. E, Normalized peak Ca²⁺ fluxes (data of D). F, ECC gain for the baseline model (solid line) and models in which dyad size and number of channels per dyad is increased 2-fold (gray line) and 3-fold (dashed line). Adapted.

Figure 6 demonstrates the ability of the Greenstein–Winslow model to reconstruct APs and cytosolic Ca²⁺ transients of normal canine midmyocardial ventricular myocytes. In Figure 6A, a normal 1-Hz steady-state AP is shown (solid black line). AP properties are similar to those measured in experiments,³⁹ with APD of ≈300 ms. Also shown in Figure 6A are the kinetics of CDI (dashed line) and VDI (dotted line), demonstrating that this model conforms with experimental findings that indicate LCC CDI is strong, whereas VDI is relatively slow and incomplete.^46,47 Figure 6B shows cytosolic (black line) and mean subspace (gray line) Ca²⁺ transients. The cytosolic Ca²⁺ transient peaks at ≈0.75 μmol/L and lasts longer than the duration of the AP whereas Ca²⁺ in the subspace reaches ≈11 μmol/L on average, and equilibrates to near-cytosolic levels rapidly, within ≈100 to 150 ms. The simulations of Figures 4 through 6 demonstrate that this multiscale model faithfully reproduces a variety of experimental observations that span from single-channel to whole myocyte.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

AP and Ca²⁺ transients at 1-Hz steady-state pacing.³⁴ A, Membrane potential (solid line, left axis) and probability that CDI and VDI (dashed and dotted lines, respectively, right axis) has not occurred as a function of time under normal conditions. B, Cytosolic (black line, left axis) and mean subspace (gray line, right axis) Ca²⁺ concentrations corresponding to the AP simulated in A. Adapted.

Additional Domains of Localized Ca²⁺ Signaling

There have been a number of recent models in which Ca²⁺ signaling in localized compartments other than the dyad have been investigated. Experimental measurements of NCX current¹⁸ have provided evidence that this transport mechanism, and possibly the sarcolemmal Ca²⁺-ATPase, may be driven by subsarcolemmal Ca²⁺, which is elevated relative to Ca²⁺ concentration in the bulk cytosol as a result of Ca²⁺ diffusion limitations. The rabbit ventricular myocyte model of Shannon et al⁵⁴ was the first to introduce a subsarcolemmal Ca²⁺ compartment, and a number of recent models now include a subsarcolemmal Ca²⁺ compartment.^35,61,92,93 Elevated subsarcolemmal Ca²⁺ may have a particularly important influence on NCX. It has been suggested that elevated subsarcolemmal Ca²⁺ during the early phase of the AP may modulate NCX reversal potential, reducing entry of Ca²⁺.⁹⁴ More needs to be done to quantify the functional consequences of sensing subsarcolemmal Ca²⁺ and NCX function. This is particularly important given the role of NCX in shaping AP properties and its possible role in the generation of delayed afterdepolarizations. Very interestingly, the Pasek et al⁹⁵ model of the guinea pig ventricular myocyte includes a diffusive t-tubule system with a heterogeneous distribution of ion channels between tubular and surface membranes, thus introducing a new extracellular compartment. This model predicted that Ca²⁺ depletion and K⁺ accumulation in the t-tubule during APs leads to a reduction in SR Ca²⁺ load, which would impact ECC.

There is now substantial electrophysiological evidence that extradyadic processes also influence CICR. Evidence suggests that Ca²⁺ entry via reverse-mode NCX activity may either trigger or modulate CICR^96,97 and that CICR modulates NCX function.⁹⁸ NCX has also been identified in the t-tubules at a high level. However, the precise location of NCX within⁹⁹ and/or near¹⁰⁰ the dyad remains uncertain. Lines et al¹⁰¹ used a model of the cardiac dyad which accounted for details of both Ca²⁺ and Na⁺ diffusion to demonstrate that if fast Na⁺ channels are in the dyadic cleft then NCX (driven by the local elevation in [Na⁺]) can influence the timing of local CICR. Recent evidence also suggests NCX may form a functional complex with NKA and its accessory protein phospholemman.¹⁰² Furthermore, one study suggests that NCX and NKA constitute a physically interacting molecular complex with ankyrin-B and an inositol 1,4,5-phosphate operated Ca²⁺-releasing channel (InsP₃R) in the t-tubules that is distinct from the LCC-RyR complex.¹⁰⁰ This hypothetical “ankyrin-B complex” may play an important role in the regulation of local Na⁺ and Ca²⁺ concentrations and possibly modulate CICR by regulating the boundary conditions on dyadic Ca²⁺. Again, it is likely that these regulatory interactions are mediated by small numbers of Ca²⁺ ions. This may be a common theme in cases where proteins interact over nanometer length scales.

Recent work also demonstrates that mitochondria are positioned near dyads.^103,104 Electron-dense structures that link the mitochondrial outer membrane with JSR, NSR, and t-tubules have been observed. These structures appear to be preferentially located, with one population forming attachments to NSR, and another set of bridges being near the dyadic cleft. Attachments to t-tubules were observed, but were rare. Ca²⁺ is taken up by mitochondria via the mitochondrial Ca²⁺-uniporter.¹⁰⁵ Therefore, spatial colocalization of mitochondria near the dyadic cleft may give rise to an additional Ca²⁺-signaling microdomain involved in the regulation of mitochondrial ATP production, and which modulates dyadic Ca²⁺ concentrations and Ca²⁺ diffusion between adjacent RyR couplons.¹⁰⁶

Regulation of CICR by Cell Signaling Pathways

Cardiac proteins involved in CICR and ECC are primary targets for regulation via cell signaling pathways and models have played a key role in understanding this regulation (see accompanying article in this series by Yang and Saucerman). The Greenstein–Winslow³⁴ model was extended to study effects of the β-adrenergic signaling pathway by the incorporation of protein kinase (PK)A-mediated functional changes in target proteins¹⁰⁷ and could explain PKA-mediated changes in AP shape, dissected the specific effects of PKA-mediated phosphorylation of LCCs and RyRs on the voltage-dependent gain of CICR, and predicted a mechanism by which increasing levels of LCC phosphorylation could lead to increased frequency of EADs.⁸⁷ Saucerman et al¹⁰⁸ developed a differential-algebraic model of the β-adrenergic signaling pathway and incorporated it into the rabbit myocyte model of Puglisi et al¹⁰⁹ and used this model to quantitatively understand and predict effects of specific molecular perturbations on cAMP dynamics and contractility via changes in function of target proteins resulting from altered phosphorylation dynamics. Recently, the model of Saucerman et al¹⁰⁸ was incorporated into the guinea pig ventricular model of Faber and Rudy¹¹⁰ in a study of the role of β-adrenergic agonists and antagonists in long-QT syndrome.¹¹¹

An important Ca²⁺ cycling regulatory mechanism of recent interest is the signaling pathway involving Ca²⁺/CaM-dependent protein kinase (CaMK)II, whose target proteins include LCCs, RyRs, phospholamban, and the SERCA pump, as well as Na⁺ and K⁺ channels.¹¹² Evidence suggests that CaMKII is directly associated with its target proteins^113,114 (implying that CaMKII signaling occurs locally in the dyadic junction), and that CaMKII activity is elevated in heart failure.¹¹⁵ The first model of the cardiac myocyte to integrate the role of CaMKII was presented by Hund and Rudy,⁴³ which predicted that CaMKII plays an important role in the rate dependent properties of the cytosolic Ca²⁺ transient, but not APD. Saucerman and Bers⁹³ incorporated models of CaM, CaMKII, and calcineurin (CaN) into the Shannon et al⁵⁴ model to better understand the functional consequences of the different affinities of CaM for CaMKII and CaN during APs. The model predicted that in the cardiac dyad, Ca²⁺ levels led to a high degree of CaM activity that results in frequency-dependent CaMKII activation and constitutive CaN activation, whereas the lower Ca²⁺ levels in the cytosol only minimally activate CaM, which allows for gradual CaN activation but no significant activation of CaMKII. Grandi et al¹¹⁶ developed a model of CaMKII overexpression in the rabbit ventricular myocyte, including the role of CaMKII phosphorylation of fast Na⁺ channels, LCCs, and the transient outward K⁺ current (I_to1), which revealed a net effect of APD reduction with CaMKII. Recently, Hashambhoy et al³⁶ described dynamic CaMKII phosphorylation of LCCs in an extension of the Greenstein–Winslow³⁴ model. In this model, it is assumed that a single CaMKII holoenzyme is tethered to each LCC, each CaMKII monomer can transition among a variety of activity states (see figure 2 of the article by Hashambhoy et al³⁶), and CaMKII monomers can catalyze phosphorylation of individual LCCs. This model demonstrated that CaMKII-dependent shifts of LCC gating patterns into high-activity gating modes may be the underlying mechanism of a variety of experimentally observed phenomena associated with I_CaL facilitation. Hashambhoy et al⁵⁵ further expanded this model to include CaMKII-dependent regulation of RyRs and demonstrated that under physiological conditions, CaMKII phosphorylation of LCCs ultimately has a greater effect on RyR leak flux and APD as compared with phosphorylation of RyRs (Figure 7A). APD was shown to correlate well with a CaMKII-mediated shift in modal gating of LCCs (Figure 7B).

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

A, Simulated results from a 1-Hz AP-pacing protocol under different LCC and/or RyR phosphorylation conditions. B, APD as a function of average LCC phosphorylation levels. Fully phosphorylated LCCs gate in mode 2, which exhibit long-duration openings.⁵⁵ Adapted.

Future Directions in ECC Modeling

A number of mechanistic questions remain unanswered in cardiac ECC. As noted above, controversy remains as to the mechanisms involved in termination of SR Ca²⁺ release. The high degree of technical difficulty involved in experimentally assessing RyR and CaRU properties under physiological conditions has resulted in limited and difficult-to-reproduce measures of certain key features of CICR such as RyR-gating properties (eg, coupled gating of RyRs), the degree of local JSR Ca²⁺ depletion associated with Ca²⁺ sparks (Ca²⁺ blinks), and the rate of local refilling of JSR. Models have been able to reveal and demonstrate multiple possible mechanisms of SR Ca²⁺ release termination, but the lack of consensus on a single experimentally verified mechanism has allowed a variety of mathematical model formulations of SR Ca²⁺ release to coexist. Improved resolution of both structural and functional imaging of ECC domain ultrastructure and localized Ca²⁺ dynamics will help to resolve these questions and to further validate and refine nano- and microscale models of ECC.^10,11 Additional areas of recent interest, where mechanistic models have yet to be fully developed, include regulation of ECC by mitochondrial Ca²⁺ dynamics (see above) and functional alteration of key ECC proteins arising from oxidative stress (eg, effects of mitochondria-derived reactive oxygen species) associated with heart failure. Recently, Tadross et al¹¹⁷ used a combination of elegant models and experiments to elucidate a mechanism underlying Ca²⁺ channel CDI, by which CaM in complex with an LCC, can sense and decode local and global Ca²⁺ signals. The implications of this mechanism on our understanding of dyadic Ca²⁺ dynamics and cardiac ECC have yet to be studied with models. Finally, the mathematical techniques used in developing models that span multiple biological scales are not adequate for all applications and themselves are an ongoing area of theoretical study (see elsewhere³⁷ for a recent review of model reduction methods). Mechanistic model reproduction of phenomena such as Ca²⁺ alternans and delayed afterdepolarizations depend on the ability of a Ca²⁺ spark to trigger neighboring release sites and form intracellular Ca²⁺ waves.⁴⁸ Although model reduction techniques have captured important features of local control of ECC (as described above), the lack of low-dimensional models that can capture Ca²⁺ wave mechanisms has made it difficult to extend such models to the tissue and whole-heart level to study how these events trigger large-scale arrhythmias. The development of models that can accomplish this task in a tractable formulation will greatly improve our ability to understand the mechanistic links between disease-related changes in ECC and arrhythmia.

Sources of Funding

Supported by NIH grants HL081427 and HL087345.

Disclosures

None.