Integrative Systems Models of Cardiac Excitation–Contraction Coupling

Excitation–contraction coupling in the cardiac myocyte is mediated by a number of highly integrated mechanisms of intracellular Ca2+ transport. The complexity and integrative nature of heart cell electrophysiology and Ca2+ cycling has led to an evolution of computational models that have played a crucial role in shaping our understanding of heart function. An important emerging theme in systems biology is that the detailed nature of local signaling events, such as those that occur in the cardiac dyad, have important consequences at higher biological scales. Multiscale modeling techniques have revealed many mechanistic links between microscale events, such as Ca2+ binding to a channel protein, and macroscale phenomena, such as excitation–contraction coupling gain. Here, we review experimentally based multiscale computational models of excitation–contraction coupling and the insights that have been gained through their application.

I ntracellular calcium (Ca 2ϩ ) concentration plays a number of important regulatory roles in the cardiac myocyte. These include excitation-contraction coupling (ECC), 1 intracellular signaling, 2 regulation of mitochondrial function, 3 gene expression, 4 and cell death. 2,5 This review focuses on the role of Ca 2ϩ in ECC in cardiac ventricular myocytes. In particular, we review experimentally based computational models of ECC, with a specific focus on mechanisms of Ca 2ϩ -induced Ca 2ϩ release (CICR) and the insights that have been gained through application of these models.
The nature of ECC is linked closely to both the microanatomical structure of the cell, and to the arrangement of contractile proteins within the cell. The basic unit of contrac-tion in the cardiac myocyte is the sarcomere. Individual sarcomeres are bounded on both ends by the t-tubular system. 1,6,7 The t-tubules are cylindrical invaginations of the sarcolemma that extend deep into the cell and approach an organelle known as the sarcoplasmic reticulum (SR). The SR is a luminal organelle located throughout the interior of the cell, and it is involved in the uptake, sequestration and release of Ca 2ϩ . SR can be divided into 3 main components known as junctional SR (JSR), corbular SR, and network SR (NSR). The JSR is that portion of the SR most closely approximating (within 12 to 15 nm) 8 the t-tubules. The close proximity of these 2 structures forms a microdomain known as the dyad. Ca 2ϩ -sensitive Ca 2ϩ -release channels known as ryanodine receptors (RyRs) are preferentially located in the dyadic region of the JSR membrane. In addition, sarcolemmal voltage-gated L-type Ca 2ϩ channels (LCCs) are preferentially located within the dyadic region of the t-tubules, where they are in close opposition to the RyRs. Early data from freeze fracture electron micrographs suggested that RyRs are arranged in a somewhat regular lattice. 9 The assumption that RyRs are tightly packed in dyadic junctions led to estimates of 100 to 300 RyRs per dyad. 8 However, more recent imaging studies have demonstrated that dyadic clefts are small, where junctions contain an average of 14 RyRs, with some clusters incompletely filled with RyRs, and with a large fraction of clusters closely spaced within 20 to 50 nm, suggesting that smaller clusters may act together to function as a single site of CICR. 10,11 During the initial depolarization stages of the action potential (AP), LCCs open allowing for the entry of Ca 2ϩ . This inward flux of Ca 2ϩ contributes to maintenance of the plateau phase of the AP. As Ca 2ϩ concentration in the dyad increases, Ca 2ϩ binds to the RyRs, increasing their open probability. Opening of RyRs leads to Ca 2ϩ release from the JSR. The amount of Ca 2ϩ released from the JSR is significantly more than the amount of trigger Ca 2ϩ entering via LCCs. 12 The ratio of release to trigger flux of Ca 2ϩ is known as ECC gain. The Ca 2ϩ release flux is a smooth, continuous function of trigger influx, a behavior observed originally by Fabiato 13 and known as graded Ca 2ϩ release. The resulting rapid, large increase in dyad Ca 2ϩ concentration ([Ca 2ϩ ] d ) leads to Ca 2ϩ diffusion from the dyad into the cytosol. As cytosolic Ca 2ϩ concentration ([Ca 2ϩ ] c ) increases, Ca 2ϩ binding to troponin C in the myofilaments initiates contraction of the myocyte.
Several mechanisms exist to remove trigger Ca 2ϩ from the myocyte, and for reuptake of released Ca 2ϩ back into the NSR. During diastole, the Na ϩ /Ca 2ϩ exchanger (NCX) is believed to import 3 Na ϩ ions for every Ca 2ϩ ion extruded, yielding a net inward charge movement. 14 (See elsewhere 15 for recently proposed alternative transport modes.) It is the principal means by which Ca 2ϩ is extruded from the myocyte. It can also function in reverse mode during the AP, in which case it extrudes Na ϩ and imports Ca 2ϩ . 16 The duration of reverse mode NCX, and hence its effect on AP shape, has been shown to vary with species, heart failure, and other factors. 17,18 The sarcolemmal Ca 2ϩ -ATPase also transports Ca 2ϩ out of the cell, but flux through this pathway is small (2% to 3%) relative to other pathways. 19 The major mechanism by which released Ca 2ϩ is transported back into the SR is the SR Ca 2ϩ -ATPase (SERCA pump). Phospholamban is normally bound to the SERCA pump, inhibiting its activity. Phosphorylation of phospholamban functionally upregulates SERCA by increasing its Ca 2ϩ affinity. 1,20 Density of the SERCA pump is highest in NSR. JSR refilling occurs via Ca 2ϩ diffusion from the NSR. As each of these mechanisms restores [Ca 2ϩ ] c to diastolic values, Ca 2ϩ dissociates from troponin C, terminating the contraction phase of the cardiac cycle.
The elucidation of CICR mechanisms has become possible with the development of experimental techniques for simul-taneous measurement of L-type Ca 2ϩ current (I CaL ) and Ca 2ϩ transients, and detection of local Ca 2ϩ release events known as Ca 2ϩ sparks. 21,22 As indicated by Soeller and Cannell, 23 the tight regulation of SR Ca 2ϩ release by triggering LCC current became evident with the following observations: (1) ECC gain is Ϸ10 in response to a voltage clamp to 0 mV 24,25 ; (2) the magnitude of SR release is a smoothly graded function of I CaL amplitude 25,26 ; (3) the bell-shaped voltage dependence of SR release is similar to that for I CaL , but with its peak shifted in the hyperpolarizing direction by Ϸ10 mV. 25,26 These experimental observations gave rise to and later verified the local control theory of ECC. The mechanism of local control predicts that within discrete regions of the dyadic cleft, the opening of an individual LCC will trigger Ca 2ϩ release from a small cluster of RyRs that reside in close proximity directly across the dyad. Tight regulation of CICR is achieved because LCCs and RyRs are regulated by [Ca 2ϩ ] d rather than by [Ca 2ϩ ] c . The amplitude and profile of the whole-cell intracellular Ca 2ϩ transient results from the recruitment and temporal summation of the elementary Ca 2ϩ spark release events. The recruitment of Ca 2ϩ sparks via this mechanism allows for graded control of SR Ca 2ϩ release. 27 Many advances in integrative modeling of the ventricular myocyte have resulted from quantitative modeling of these components which underlie CICR and ECC in cardiac myocytes. In this review, we will present a range of models that have been developed to understand properties of CICR and ECC in the cardiac ventricular myocyte. We have chosen to focus on development of major new functional components of CICR and ECC, their integration into whole-cell models, and the insights that have been obtained from these advances in modeling capabilities that span multiple biological scales.

Common Pool Models of CICR
The landmark DiFrancesco-Noble model of the cardiac Purkinje fiber 28 was developed from the preceding Beeler-Reuter model 29 (first computational model of the cardiac ventricular myocyte) and that of McAllister et al. 30 This model was ground breaking for many reasons and was the first model to introduce descriptions of the transport of Ca 2ϩ into the NSR by SERCA, extrusion of Ca 2ϩ from the cytosol by NCX, diffusion of Ca 2ϩ 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 2ϩ 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 2ϩ release first lead to a series of "common pool" models, which, in fact, fail to capture these fundamental properties. As defined by Stern, 27 common pool models are those in which LCC trigger Ca 2ϩ enters the same Ca 2ϩ pool into which JSR Ca 2ϩ is released, and which controls JSR Ca 2ϩ release. The rapid increase of Ca 2ϩ in this pool leads to nonphysiological regenerative, all-or-none rather than graded Ca 2ϩ release, a situation where whole-cell SR Ca 2ϩ release escapes control by I CaL once it is initiated. 26 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 2ϩ release. 27 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 31 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 32 based on single-LCC recordings in which Ca 2ϩ 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 2ϩ 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][34][35] as well as other modes such as phosphorylation. 36 The Jafri et al 31 model was also first to implement a "restricted subspace," a single compartment representing the total volume of all dyads into which all Ca 2ϩ flux through RyRs and LCCs is directed. This model made a number of predictions regarding mechanisms of CICR and Ca 2ϩ cycling (reviewed elsewhere 37 ). 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 2ϩ load dynamics. This Ca 2ϩ subsystem was modified for canine in the model of Winslow et al 38 from which a model of heart failure was developed based on data obtained from myocytes isolated from dogs with tachycardia pacing-induced heart failure. 39 The model predicted that in heart failure, downregulation of the SERCA pump and upregulation of NCX decreased JSR Ca 2ϩ load, Ca 2ϩ 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 2ϩ 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  31 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. 73 Reprinted with permission from Elsevier. degeneration of electric propagation into a cardiac arrhythmia, whereas the control myocyte model yielded normal coordinated ventricular activation during pacing. 42 Despite the fact that these common pool models were not able to reproduce graded SR Ca 2ϩ 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 2ϩ 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 2ϩ release. The critical relationship between SR Ca 2ϩ release and APD, 38,40,41 and the recognition that the feedback of SR Ca 2ϩ release on LCCs occurs locally in discrete Ca 2ϩ spark events 21 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][46][47] Strong CDI in combination with weak/incomplete VDI causes instability in common pool models 34 (which were therefore typically formulated with strong VDI and weak CDI) because all-or-none SR Ca 2ϩ 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 34 ). 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 2؉ Release
In 1992, Stern 27 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 2ϩ released from the SR from the trigger Ca 2ϩ 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 2ϩ ions permeating an LCC privileged access to the Ca 2ϩ binding sites on nearby RyRs. Importantly, the assumed structural separation of individual clefts 8 means that released Ca 2ϩ 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 2ϩ signal. In this model, graded release is achieved by statistical recruitment of elementary SR Ca 2ϩ 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 2ϩ distribution within the dyadic junction, some of which focused on mechanisms of Ca 2ϩ 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 2ϩ release units.
Because of the small size of the dyad, it became clear that understanding the spatial and temporal profile of [Ca 2ϩ ] d would be an important step in elucidating the mechanistic details of local CICR. Langer and Peskoff 49 developed one of the first detailed models to numerically simulate the time course of dyadic Ca 2ϩ . The model predicted that a single typical LCC opening led to Ϸ1 mmol/L [Ca 2ϩ ] d near the mouth of the LCC. The presence of sarcolemmal buffering sites slowed the decay of the local Ca 2ϩ 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 50 formulated a dyad model that introduced the electrostatic effects on ion movement that predicted that [Ca 2ϩ ] d actually changes more rapidly in response to LCC gating. In this model, [Ca 2ϩ ] d was approximately proportional to unitary LCC currents over a wide range of current amplitudes. The model yielded predictions of [Ca 2ϩ ] d profiles during JSR release, which peaked as high as 300 mol/L near the center of the dyad. The decay in dyadic Ca 2ϩ 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 2ϩ 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 2ϩ sparks are the elementary events underlying SR Ca 2ϩ release (for a review see 51 ), mathematical models of Ca 2ϩ sparks were developed. The summation of Ca 2ϩ sparks could explain the characteristic properties of CICR, ECC gain, and whole-cell Ca 2ϩ transients; therefore, the analysis of models of single Ca 2ϩ sparks and ensembles of Ca 2ϩ sparks provided a wealth of information. Local control of CICR requires that SR release be stable, which suggests that Ca 2ϩ sparks, which are locally regenerative, must be terminated rapidly by a robust mechanism. Although the mechanisms of Ca 2ϩ spark termination still remain unclear, a number of models have helped elucidate how possible mechanisms may contribute to termination of SR Ca 2ϩ release.
Stern 27 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 2ϩ 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) 52 or taking on the complex properties of adaptation. 51 (2) Depletion of locally releasable JSR Ca 2ϩ will occur as the result of RyR activity. The emptying of the cisternal Ca 2ϩ store reduces the Ca 2ϩ gradient across the RyRs and hence the Ca 2ϩ 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 2ϩ diffuses rapidly out of the dyad, the simultaneous closure of all RyRs would lead to a rapid reduction in [Ca 2ϩ ] 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 2ϩ spark, and it does not provide for a mechanism of refractoriness. The stereotypical nature of Ca 2ϩ 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. 23 An additional mechanism proposed recently is RyR regulation by luminal SR Ca 2ϩ . RyR sensitivity to [Ca 2ϩ ] d can be altered by the level of JSR luminal [Ca 2ϩ ] via the interaction of the Ca 2ϩ buffer calsequestrin with RyR accessory proteins triadin and junction within the JSR. 53 The decrease in RyR open probability with JSR depletion likely contributes to termination of Ca 2ϩ release.
In 1999, Stern et al 33 developed a stochastic model of cardiac ECC. The model consisted of a cardiac dyad in which the Ca 2ϩ concentration profile was simulated by solving equations governing diffusion and Ca 2ϩ 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 2ϩ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 33 ). The work of Stern et al 33 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 2ϩ release remains controversial (see below). A number of RyR models 57,58 were developed to better understand the observation of RyR adaptation, 59 whereby an RyR can reopen in response successive increases in Ca 2ϩ concentration, indicating they have not entered a conventional absorbing CDI state. Rice 33 which introduced an important modulatory role for SR luminal Ca 2ϩ 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 54 and human, 61 the role of RyR phosphorylation on Ca 2ϩ dynamics and the AP, 55 and the role of spatial coupling between dyadic release sites in Ca 2ϩ spark-induced sparks and Ca 2ϩ alternans. 48 In 2002, Sobie et al 62 developed a Ca 2ϩ spark model to further investigate mechanisms underlying termination of SR Ca 2ϩ 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 2ϩ increases the Ca 2ϩ -sensitivity of RyR activation 63 ; and (2) RyR coupling cooperatively influences their gating properties within a cluster. 64 These mechanisms of RyR gating yield a "sticky cluster" model which reproduces fundamental features of Ca 2ϩ 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 33 ) 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 62 model may not be realistic as experiments have demonstrated that 40% to 75% of releasable Ca 2ϩ 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 2ϩ rather than an RyR inactivation mechanism.
Although many controversies as to the details of CICR remain to be resolved, the models of dyadic Ca 2ϩ dynamics, Ca 2ϩ 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 2ϩ 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 34 (Greenstein-Winslow model). This model is described in detail in the following section.

Multiscale Modeling of CICR
Quantitative understanding of CICR requires that Ca 2ϩ 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][72][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 2ϩ signaling in the dyad that simulates motion of Ca 2ϩ ions, binding of ions to RyR binding sites, and Ca 2ϩ 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 2ϩ 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 2ϩ 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 2ϩ ions, as influenced by thermal energy and electrostatic potentials.
The small size of the dyadic cleft (Ϸ4ϫ10 5 nm 3 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 2ϩ 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 2ϩ concentration in the dyad ranges from 100 to 1000 mol/ L, 49,50 which corresponds to as few as 20 to 200 free Ca 2ϩ ions. 73 Therefore, application of laws of mass-action to analysis of signaling events mediated by this relatively small number of ions may not be justified. 77 Rather, the stochastic motion of Ca 2ϩ ions in the dyad may impart a degree of signaling noise between LCCs and RyRs, thereby influencing properties of CICR. 73 In addition, the physical location and dimensions of dyad proteins may have considerable influence on motion of Ca 2ϩ ions and hence, CICR.
Development of the nanoscale model is presented in detail elsewhere. 73 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 2ϩ ions into the dyad via open LCCs and RyRs was modeled as a Poisson process with rate governed by channel permeabilities measured experimentally. 79 CDI of an LCC required that the 2 carboxylterminal Ca 2ϩ binding sites of the associated CaM were occupied. Opening of an RyR required that Ca 2ϩ was bound to all 4 activation sites. RyRs were assumed to inactivate when Ca 2ϩ was bound to at least 1 of 4 distinct inactivation sites. Sarcolemmal anionic sites acted as Ca 2ϩ buffers. 49,50 The joint positions of Ca 2ϩ ions present in the dyad were modeled as a Brownian motion in a potential field, describing both interactions of Ca 2ϩ ions with other Ca 2ϩ ions as well as electrostatic potentials based on the Fokker-Planck equation (FPE). 80 We generated sample paths of Ca 2ϩ 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. 81 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 2ϩ ions within a 50 nm radius of the LCC (see figure 8 in the article by Tanskanen et al 73 ), with a minimum of 0 and a maximum of 7 Ca 2ϩ ions, corresponding to 98.7 mol/L [Ca 2ϩ ] 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 2ϩ ions. The number of Ca 2ϩ ions peaked at 49 during Ca 2ϩ release, which corresponds to a concentration of Ϸ0.7 mmol/L Ca 2ϩ , similar to the concentration of Ca 2ϩ found in the JSR. 82 With a coefficient of variation of 64%, RyR gating produces a high degree of variability in the number of dyadic Ca 2ϩ ions, showing that the LCC-RyR signaling involved in CICR at the dyad level is a noisy process mediated by tens of Ca 2ϩ ions. Figure 2A and 2B demonstrates fundamental features of a representative Ca 2ϩ 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 2ϩ 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 2ϩ 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 2ϩ signal in the dyad (Figure 2A) is closely correlated to the number of open RyRs ( Figure 2B). Figure 2C shows the average dyadic Ca 2ϩ signal, calculated by averaging the number of Ca 2ϩ ions at each point in time across 400 independent dyads. The duration of the average local Ca 2ϩ transient shown in Figure  2C is Ϸ20 ms (at half-maximal amplitude), similar to that measured for Ca 2ϩ spikes in myocytes. 21,24,52 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 73 ). 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 2ϩ 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 2ϩ -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 2ϩ 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 2ϩ 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 2ϩ 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 34 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 38 to include a population of dyadic Ca 2ϩ 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 2ϩ 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 2ϩ 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 2ϩ 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. 34 Figure 3A shows a schematic of the CaRU model which is intended to mimic the properties of Ca 2ϩ 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. 83 The division of the CaRU into 4 subunits allows for the possibility that an LCC may trigger Ca 2ϩ release in adjacent subspaces (ie, RyR recruitment). On activation of RyRs, subspace [Ca 2ϩ ] will become elevated. This Ca 2ϩ 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 2ϩ 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 measure- ments in canine myocytes, 84 which corresponded to 50 000 LCCs (12 500 CaRUs) and agrees with experimental estimates of LCC density. 79 LCC and RyR gating were based on previous models, 31,85 and the channels (ClCh) that carry the Ca 2ϩ -dependent transient outward chloride (Cl -) current (I to2 ) are included within the CaRU. 86 A detailed description of the local control simulation algorithm is given in Appendix II of a previously published article. 34 Figure 4B shows the corresponding subspace [Ca 2ϩ ], which reaches a peak value of Ϸ40 mol/L. Total Ca 2ϩ influx through the set of 4 LCCs (gray line) and the net SR Ca 2ϩ 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 2ϩ 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 2ϩ concentration in the 4 subspace compartments of the CaRU is intended to represent a Ca 2ϩ 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 Previous experimental studies 25,88 and mathematical models 27,33 have shown significant differences between the voltage dependence of I CaL and RyR Ca 2ϩ release flux even though SR Ca 2ϩ release is exclusively controlled by LCC Ca 2ϩ 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 2ϩ flux. In Figure 5B, the peak fluxes of Figure 5A have been normalized based on their respective maxima. Maximal LCC Ca 2ϩ influx occurs at 10 mV, whereas maximal RyR Ca 2ϩ 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 volt-age 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. Figure 6 demonstrates the ability of the Greenstein-Winslow model to reconstruct APs and cytosolic Ca 2ϩ 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, 39 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 2ϩ transients. The cytosolic Ca 2ϩ transient peaks at Ϸ0.75 mol/L and lasts longer than the duration of the AP whereas Ca 2ϩ 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 singlechannel to whole myocyte.

The Coupled LCC-RyR Local Control Model
A limitation of the Greenstein-Winslow 34 local control model is that it is far more computationally demanding than deterministic models. To address this, we recently developed an analytic "first-principles" approach for deriving simplified mechanistic models of CICR by applying a carefully chosen set of approximations that allow for the ensemble behavior of CaRUs to be represented by a low dimensional system of ODEs, eliminating the need for computationally expensive stochastic simulations. This model, termed the coupled LCC-RyR model, 72 retains the biophysically based description of local control of CICR, and captures its key properties including graded SR Ca 2ϩ release and voltage dependent ECC gain and is the basis of a deterministic local control model of the cardiac myocyte. 71 Following the methods described by Hinch et al, 72 LCC and RyR models composed of the minimum number of states necessary to capture fundamental dynamics of channel gating were designed and constrained based on experimental data sets. A simplified CaRU model consisting of only a single LCC and a single RyR was formed. The relatively fast diffusion (ie, short diffusion distance) of Ca 2ϩ from dyad to cytosol (as compared to the time scale of LCC and RyR gating) 62 allows for the application of a steady-state approximation for dyadic Ca 2ϩ diffusion. This allows the resulting CaRU to be analytically simplified into a single Markov state model in which each state represents a unique configuration of channel states for the combined set of LCCs and RyRs in the CaRU. In the resulting formulation, the dynamics of a population of CaRUs can be described by a deterministic set of coupled ODEs. Other elegant methods have been derived to improve computational efficiency of ECC simulations as well, 37,89,90 but have not yet been implemented in an integrated myocyte model using biophysically detailed channel models. Figure 5D through 5F demonstrates that the deterministic coupled LCC-RyR model reproduces the fundamental quantitative properties of CICR in a similar manner to its stochastic predecessor (in Figure 5A through 5C). 34 The effect of dyad size on ECC was analyzed using models with 2-fold and 3-fold larger dyads ( Figure 5F, gray and dashed lines, respectively). In each case, the number of dyads per cell was reduced such that the total number of LCCs and RyRs per cell was conserved, and the RyR/LCC stoichiometry was unchanged. At potentials more negative than Ϫ10 mV, the increase in CaRU size yields a substantial increase in gain such that the ECC gain curve acquires the steep slope similar to those observed experimentally. 24,25 This arises as a result of a model-dependent increase in the functional RyR/LCC ratio because it is likely that a single LCC opening will trigger regenerative release from most RyRs in a CaRU. The coupled LCC-RyR local control model reconstructs APs and Ca 2ϩ transients of canine midmyocardial ventricular cells (see figure 6 of the article by Greenstein et al 71 for details) at a variety of pacing frequencies, and similar to the Greenstein-Winslow model, has an I CaL which exhibits strong CDI and weak VDI. The ability of this model to reproduce all of the core features of the stochastic local control model in an efficient manner has made it a candidate for use in large scale tissue simulations. 91

Additional Domains of Localized Ca 2؉ Signaling
There have been a number of recent models in which Ca 2ϩ signaling in localized compartments other than the dyad have been investigated. Experimental measurements of NCX current 18 have provided evidence that this transport mechanism, and possibly the sarcolemmal Ca 2ϩ -ATPase, may be driven by subsarcolemmal Ca 2ϩ , which is elevated relative to Ca 2ϩ concentration in the bulk cytosol as a result of Ca 2ϩ diffusion limitations. The rabbit ventricular myocyte model of Shannon et al 54 was the first to introduce a subsarcolemmal Ca 2ϩ compartment, and a number of recent models now include a subsarcolemmal Ca 2ϩ compartment. 35,61,92,93 Elevated subsarcolemmal Ca 2ϩ may have a particularly important influence on NCX. It has been suggested that elevated subsarcolemmal Ca 2ϩ during the early phase of the AP may modulate NCX reversal potential, reducing entry of Ca 2ϩ . 94 More needs to be done to quantify the functional consequences of sensing subsarcolemmal Ca 2ϩ 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 95 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 2ϩ depletion and K ϩ accumulation in the t-tubule during APs leads to a reduction in SR Ca 2ϩ load, which would impact ECC.
There is now substantial electrophysiological evidence that extradyadic processes also influence CICR. Evidence suggests that Ca 2ϩ entry via reverse-mode NCX activity may either trigger or modulate CICR 96,97 and that CICR modulates NCX function. 98 NCX has also been identified in the t-tubules at a high level. However, the precise location of NCX within 99 and/or near 100 the dyad remains uncertain. Lines et al 101 used a model of the cardiac dyad which accounted for details of both Ca 2ϩ 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. 102 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 2ϩ -releasing channel (InsP 3 R) in the t-tubules that is distinct from the LCC-RyR complex. 100 This hypothetical "ankyrin-B complex" may play an important role in the regulation of local Na ϩ and Ca 2ϩ concentrations and possibly modulate CICR by regulating the boundary conditions on dyadic Ca 2ϩ . Again, it is likely that these regulatory interactions are mediated by small numbers of Ca 2ϩ 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 2ϩ is taken up by mitochondria via the mitochondrial Ca 2ϩ -uniporter. 105 Therefore, spatial colocalization of mitochondria near the dyadic cleft may give rise to an additional Ca 2ϩ -signaling microdomain involved in the regulation of mitochondrial ATP production, and which modulates dyadic Ca 2ϩ concentrations and Ca 2ϩ diffusion between adjacent RyR couplons. 106

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 34 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 107 and could explain PKAmediated 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. 87 111 An important Ca 2ϩ cycling regulatory mechanism of recent interest is the signaling pathway involving Ca 2ϩ /CaMdependent protein kinase (CaMK)II, whose target proteins include LCCs, RyRs, phospholamban, and the SERCA pump, as well as Na ϩ and K ϩ channels. 112 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. 115 The first model of the cardiac myocyte to integrate the role of CaMKII was presented by Hund and Rudy, 43 which predicted that CaMKII plays an important role in the rate dependent properties of the cytosolic Ca 2ϩ transient, but not APD. Saucerman and Bers 93 Figure 7A). APD was shown to correlate well with a CaMKII-mediated shift in modal gating of LCCs ( Figure 7B).

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 2ϩ 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 2ϩ depletion associated with Ca 2ϩ sparks (Ca 2ϩ blinks), and the rate of local refilling of JSR. Models have been able to reveal and demonstrate multiple possible mechanisms of SR Ca 2ϩ release termination, but the lack of consensus on a single experimentally verified mechanism has allowed a variety of mathematical model formulations of SR Ca 2ϩ release to coexist. Improved resolution of both structural and functional imaging of ECC domain ultrastructure and localized Ca 2ϩ 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 2ϩ 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 117 used a combination of elegant models and experiments to elucidate a mechanism underlying Ca 2ϩ channel CDI, by which CaM in complex with an LCC, can sense and decode local and global Ca 2ϩ signals. The implications of this mechanism on our understanding of dyadic Ca 2ϩ 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 37 for a recent review of model reduction methods). Mechanistic model reproduction of phenomena such as Ca 2ϩ alternans and delayed afterdepolarizations depend on the ability of a Ca 2ϩ spark to trigger neighboring release sites and form intracellular Ca 2ϩ waves. 48 Although model reduction techniques have captured important features of local control of ECC (as described above), the lack of lowdimensional models that can capture Ca 2ϩ 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 diseaserelated changes in ECC and arrhythmia.

Conclusion
Integrative modeling of cardiac ECC and myocyte physiology has played a critical role in revealing mechanistic insights at and across a range of biological scales. At the smallest scale, models have shed light on how Ca 2ϩ ions and proteins in the cardiac junction interact during LCC and RyR gating. Expanding our view, simplified models of individual release sites have allowed for detailed simulations that reveal determinants of Ca 2ϩ spark activation and termination properties. On an intermediate scale, models of independent dyadic junctions have revealed how characteristic whole-myocyte properties of ECC emerge from the ensemble behavior of individual release sites. Incorporation of such Ca 2ϩ release sites into whole-cell models elucidates the bidirectional interaction between ECC and membrane potential that determine properties of the cardiac AP. Incorporation of ECC regulatory signaling pathways into models has further refined our understanding of this integrative system (see accompanying article in this series by Yang and Saucerman). At the highest scale, such models can be used as building blocks for large-scale tissue and whole-heart models that can reveal how wave propagation in a normal heart tissue can degenerate into an arrhythmia as a result of defective ECC (see accompanying articles in this series by Weiss et al and Trayanova et al). The continuum of biological scales spanned by these models leads to the formation of powerful multiscale approaches whereby macroscopic phenotypes (eg, steep ECC gain curve) can be understood as resulting from microscopic mechanisms (eg, Ca 2ϩ funneling by RyR protein geometry). In the same way that Hodgkin and Huxley formed a mathematical framework that now defines the way in which we conceptualize ion channel function, modern multiscale models of cardiac ECC have provided a framework for understanding complex ECC mechanisms, interpreting experimental data, and making mechanistic predictions that guide experimental design.

Sources of Funding
Supported by NIH grants HL081427 and HL087345.