Copyright © 2008 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 95, Issue 4, 1689-1703, 15 August 2008
doi:10.1529/biophysj.107.125948
Biophysical Theory and Modeling
Moment Closure for Local Control Models of Calcium-Induced Calcium Release in Cardiac Myocytes
George S.B. Williams*, 1, Marco A. Huertas*, Eric A. Sobie†, M. Saleet Jafri‡ and Gregory D. Smith*, §, 
, 
* Department of Applied Science, College of William and Mary, Williamsburg, Virginia
† Department of Bioinformatics and Computational Biology, George Mason University, Manassas, Virginia
‡ Department of Pharmacology and Systems Therapeutics, Mount Sinai School of Medicine, New York, New York
§ Mathematical Biosciences Institute, The Ohio State University, Columbus, Ohio
Address reprint requests to Gregory D. Smith.
1 George S. B. Williams and Marco A. Huertas contributed equally to this work.
Abstract
In prior work, we introduced a probability density approach to modeling local control of Ca
2+-induced Ca
2+ release in cardiac myocytes, where we derived coupled advection-reaction equations for the time-dependent bivariate probability density of subsarcolemmal subspace and junctional sarcoplasmic reticulum (SR) [Ca
2+] conditioned on Ca
2+ release unit (CaRU) state. When coupled to ordinary differential equations (ODEs) for the bulk myoplasmic and network SR [Ca
2+], a realistic but minimal model of cardiac excitation-contraction coupling was produced that avoids the computationally demanding task of resolving spatial aspects of global Ca
2+ signaling, while accurately representing heterogeneous local Ca
2+ signals in a population of diadic subspaces and junctional SR depletion domains. Here we introduce a computationally efficient method for simulating such whole cell models when the dynamics of subspace [Ca
2+] are much faster than those of junctional SR [Ca
2+]. The method begins with the derivation of a system of ODEs describing the time-evolution of the moments of the univariate probability density functions for junctional SR [Ca
2+] jointly distributed with CaRU state. This open system of ODEs is then closed using an algebraic relationship that expresses the third moment of junctional SR [Ca
2+] in terms of the first and second moments. In simulated voltage-clamp protocols using 12-state CaRUs that respond to the dynamics of both subspace and junctional SR [Ca
2+], this moment-closure approach to simulating local control of excitation-contraction coupling produces high-gain Ca
2+ release that is graded with changes in membrane potential, a phenomenon not exhibited by common pool models. Benchmark simulations indicate that the moment-closure approach is nearly 10,000-times more computationally efficient than corresponding Monte Carlo simulations while leading to nearly identical results. We conclude by applying the moment-closure approach to study the restitution of Ca
2+-induced Ca
2+ release during simulated two-pulse voltage-clamp protocols.
Introduction
The key step linking electrical excitation to contraction in cardiac myocytes is Ca2+-induced Ca2+ release (CICR), in which Ca2+ current flowing across the cell membrane triggers the release of additional Ca2+ from the sarcoplasmic reticulum (SR). In ventricular cells, CICR occurs as a set of discrete microscopic events known as Ca2+ sparks 1, with each spark triggered by local, rather than cell-wide, increases in myoplasmic [Ca2+]. As a consequence of this local-control mechanism of CICR, the cellular SR Ca2+ release flux is not a function of a single quantity, such as spatially averaged intracellular [Ca2+], but instead depends on thousands of different local Ca2+ concentrations, each of which can fluctuate with stochastic openings and closings of nearby Ca2+ channels in the sarcolemmal and SR membranes. The picture is further complicated by the fact that dynamic changes in local SR [Ca2+], which are also spatially heterogeneous, are thought to influence the gating of SR Ca2+ release channels known as ryanodine receptors (RyRs).
Computational models have been developed in which SR Ca2+ release depends directly on the average myoplasmic [Ca2+] 2,3,4. These so-called common-pool models 5 display SR Ca2+ release that occurs in an all-or-none fashion, contrary to experiments showing that release is smoothly graded with changes in Ca2+ influx 6,7,8. On the other hand, several published models achieve graded Ca2+ release using nonmechanistic formulations, such as having SR Ca2+ release depend explicitly on Ca2+ currents rather than on local [Ca2+] 9,10,11,12,13.
Models of EC coupling are able to reproduce graded Ca2+ release mechanistically by simulating the stochastic gating of channels in Ca2+ release sites using Monte Carlo methods. In these approaches, one or more L-type Ca2+ channels interact with a cluster of RyRs through changes in [Ca2+] in a small diadic subspace between the sarcolemmal and SR membranes. These models also generally consider local changes in junctional SR [Ca2+], because these changes are thought to be important for Ca2+ spark termination and refractoriness 14,15,16. Realistic cellular SR Ca2+ release can be simulated by computing the stochastic triggering of sparks from hundreds to thousands of such Ca2+ release units (CaRUs) 5,15,16,17. However, Monte Carlo simulations of local control of EC coupling can be computationally demanding, making it difficult to augment these models with representations of the ionic currents responsible for action potentials, and impractical to use this approach for simulations of phenomena occurring over the course of many heartbeats.
We recently demonstrated that an alternative probability-density approach can be used to simulate graded, locally controlled SR Ca2+ release mechanistically 18. In this prior work, coupled advection-reaction equations were derived relating the time-dependent probability density of subsarcolemmal subspace and junctional SR [Ca2+] conditioned on CaRU state. By numerically solving these equations using a high-resolution finite difference scheme and coupling the resulting probability densities to ordinary differential equations (ODEs) for the bulk myoplasmic and sarcoplasmic reticulum [Ca2+], a realistic but minimal model of cardiac excitation-contraction coupling was produced. This new approach to modeling local control of EC coupling is often computationally more efficient than Monte Carlo simulation, particularly if the dynamics of subspace [Ca2+] are much faster that those of junctional SR [Ca2+], allowing the bivariate probability density functions for subspace and junctional SR [Ca2+] to be replaced with univariate densities for junctional SR [Ca2+]. However, the probability density approach can lose its computational advantage when the number of states in the CaRU model is large or the dynamics of local [Ca2+] are such that numerical stability requires a refined mesh for solving the advection-reaction equations.
We therefore aimed to develop methods for improving upon the probability-density approach, and in this study, we describe a moment-closure technique that leads to significant computational advantages. After briefly reviewing the Monte Carlo and probability density approaches to modeling local control of EC coupling in cardiac myocytes, the new methodology begins with a derivation of a system of ODEs describing the time-evolution of the moments of the univariate probability density functions for junctional SR [Ca2+] jointly distributed with CaRU state. This open system of ODEs is then closed using an algebraic relationship that expresses the third moment of junctional SR [Ca2+] in terms of the first and second moments. In this manner, the partial differential equations describing the univariate probability densities of junctional SR [Ca2+] jointly distributed with CaRU state are replaced with ODEs describing the time-evolution of the moments of these distributions. In simulated voltage-clamp protocols using 12-state CaRUs that respond to the dynamics of both subspace and junctional SR [Ca2+], this moment-closure approach to simulating local control of EC coupling produces high-gain Ca2+ release that is graded with changes in membrane potential, a phenomenon not exhibited by common pool models. Benchmark simulations indicate that this moment-closure technique for local control models of CICR in cardiac myocytes is nearly 10,000-times more computationally efficient than corresponding Monte Carlo simulations, while leading to nearly identical results. We conclude by applying the moment-closure approach to study the restitution of Ca2+-induced Ca2+ release during simulated two-pulse voltage-clamp protocols.
Model formulation
The focus of this article is a moment-closure technique to modeling local control of CICR in cardiac myocytes. The whole cell model of EC coupling that will be used to demonstrate the method closely follows our prior work in which we presented traditional Monte Carlo simulations of graded, locally controlled SR Ca2+ release to validate a novel probability density approach that represents the distribution of diadic subspace and junctional SR Ca2+ concentrations with a system of partial differential equations 18. Below we briefly review the Monte Carlo and probability density formulations, emphasizing minor adjustments that were required to implement the moment-closure technique. The Results section begins with the derivation of the moment-closure equations and follows with the validation and benchmarking of the moment-closure technique for local control models of CICR in cardiac myocytes by comparison to Monte Carlo simulation.
Monte Carlo formulation
The Monte Carlo model of local control of CICR in cardiac myocytes describes the dynamics of bulk myoplasmic [Ca2+], network SR [Ca2+], N diadic subspace Ca2+ concentrations, and N junctional SR domain Ca2+ concentrations through a system of ODEs. These are coupled to N Markov chains representing the stochastic gating of each CaRU that consists of one L-type Ca2+ channel (DHPR) and one RyR megachannel coupled through the local diadic subspace (cds) [Ca2+]. While a complete description of CICR would include stochastic gating of roughly N=10,000 CaRUs, each containing multiple L-type Ca2+ channels (1–10) 19 and RyRs (30–300) 20, Monte Carlo simulations of EC coupling focusing on local control have often used Markov models of reduced complexity 5,16,21. This level of resolution will suffice to introduce the moment-closure technique.
Concentration balance equations
The Monte Carlo model consists of N+2 ODEs representing the time-evolution of [Ca2+] in the bulk myoplasm (cmyo), network SR (cnsr), and N junctional SRs 
compartments. Consistent with Fig. 1, the concentration balance equations for these compartments are
Concentration balance equations
The Monte Carlo model consists of N+2 ODEs representing the time-evolution of [Ca2+] in the bulk myoplasm (cmyo), network SR (cnsr), and N junctional SRs 
compartments. Consistent with Fig. 1, the concentration balance equations for these compartments are
where 1≤
n≤
N and
λnsr and
λjsr are volume fractions (see the
Appendix ). The flux through the RyR megachannel associated with the
nth CaRU
is given by
where
is a stochastic variable that takes the value 1 or 0 depending on whether the
nth RyR megachannel is open or closed, and
is the associated diadic subspace concentration defined below (Eq.
(9)). Similarly, diffusion from the network SR to each junctional SR compartment is given by
The total refill flux occurring in Eq. (2) includes the contribution from each CaRU and is given by
while the total flux out of the
N diadic subspaces is given by
The remaining four fluxes that appear in Eqs.
(1) and
Fig. 1 include
(influx into the diadic subspaces via L-type Ca
2+ channels which are functions of the random variable
),
Jin (background Ca
2+ influx),
Jncx (Na
+-Ca
2+ exchange),
Jserca (SR Ca
2+-ATPases), and
Jleak (the network SR leak). The functional form of these four fluxes can be found in the
Appendix .
Diadic subspace Ca2+ concentration
Note that a concentration balance equation is not included for diadic subspace [Ca2+], because in our previous study we observed that model parameters lead to rapid equilibrium of the diadic subspace [Ca2+] with the [Ca2+] in the junctional SR and bulk myoplasm 18. Thus, in each diadic subspace we assume a [Ca2+] 
that balances the fluxes in and out of that compartment,
that is,
where 1≤
n≤
N and
In these expressions, the quantities
and
indicate whether the channel is open or closed, and
,
, and
are functions of plasma membrane voltage defined by
where the L-type Ca
2+ channel flux,
is given by Eq.
(56) (see the
Appendix ).
Twelve-state CaRU model
The RyR model used here is similar to the two-state minimal model of an RyR megachannel used in prior work 18. Consistent with several studies indicating that the gating of the RyR cluster associated with each CaRU is essentially all-or-none 5,15,16, the two-state RyR megachannel model used in Williams et al. 18 included transition rates that were nonlinear functions of diadic subspace (cds) and junctional SR (cjsr) [Ca2+], thereby allowing for Ca2+-dependent activation of RyR gating as well as spark termination facilitated by localized depletion of junctional SR [Ca2+]. Because the moment-closure approach is most easily presented when all Ca2+-mediated transitions in the CaRU model are bimolecular association reactions, the six-state RyR megachannel model used here employs sequential binding of diadic subspace Ca2+ ions to achieve highly cooperative Ca2+-dependent opening of the RyR megachannel. Similarly, an explicit junctional SR Ca2+-dependent transition is included so that depletion of luminal Ca2+ decreases the open probability of the megachannel,
Parameters were chosen (see
Table 3) so that the behavior of this minimal six-state RyR megachannel model approximated the above-mentioned two-state model.
As in prior work 18, we use a two-state model of the L-type Ca2+ channel (DHPR),
where
C and
represent closed and open states,
is the voltage-dependent activation rate
10 given by
and
is the constant deactivation rate that sets the mean open time (0.2ms) and maximum open probability (0.1) of the channel. Although this two-state DHPR model ignores voltage- and Ca
2+-dependent inactivation of L-type Ca
2+ channels, these processes do not significantly influence the triggering of CICR during the whole-cell voltage clamp protocols that are used in this article to validate the moment-closure technique.
Combining the six-state RyR megachannel model with the two-state L-type channel model yields a 12-state CaRU model that takes the form
where horizontal and vertical transitions are governed by Eqs.
(13), respectively, and the first character
indicates the state of the DHPR while the second character
refers to the state of the RyR megachannel.
Note that the 12×12 infinitesimal generator matrix (sometimes called the Q-matrix) that collects the rate constants of the CaRU model (Eq. (16)) can be written compactly in the form
where the elements of
are the Ca
2+-independent transitions (both voltage-dependent and voltage-independent with units of time
−1), and the elements of
Kds and
Kjsr are the association rate constants for the transitions mediated by diadic subspace (
cds) and junctional SR (
cjsr) [Ca
2+], respectively (with units of concentration
−1 time
−1). Although noncooperative binding of Ca
2+ is not a formal requirement for the application of the moment-closure technique, for simplicity we will assume the CaRU model is written in the form of Eq.
(17).
Univariate probability density model
The moment-closure technique begins with the equations for a univariate probability density model of local control of Ca2+-induced Ca2+ release in cardiac myocytes 18. We write ρi(cjsr, t) to denote probability density functions for the distribution of [Ca2+] in a large number of junctional SR compartments jointly distributed with CaRU state, that is,
where
i is an index over CaRU state, and the tilde in
and
indicate random quantities. For these densities to be consistent with the dynamics of the Monte Carlo model of cardiac EC coupling as
N→∞, they must satisfy a system of advection-reaction equations of the form
18,
22,
23where 1≤
i≤
M,
M=12 is the number of states in the CaRU model,
Q is the
M×
M generator matrix (Eq.
(17)), the row-vector
ρ(
cjsr,
t)=(
ρ1,
ρ2, ···,
ρM) collects the time-dependent probability densities for the junctional SR [Ca
2+] jointly distributed with CaRU state (Eq.
(18)), and [
ρQ]
i is the
ith element of the vector-matrix product
ρQ.
Note that the factor 
in Eq. (19) describes the deterministic aspect of the time-evolution of cjsr when the CaRU is in state i. That is, consistent with Eq. (3) we have
where 1≤
i≤
M and
is a function of CaRU state, the local junctional SR [Ca
2+], and the bulk myoplasmic [Ca
2+] analogous to Eqs.
(9),
where
In these expressions, the quantities
and
take values of 0 or 1 depending on whether the respective component of the CaRU model is closed or open, and
and
are functions of plasma membrane voltage defined by
where
is the total flux through the L-type Ca
2+ channels (Eq.
(58)).
Conversely, the reaction terms ([ρQ]i) on the right-hand side of Eq. (19) correspond to the stochastic aspect of the CaRU dynamics (i.e., changes in probability due to the stochastic gating of the RyR megachannel and DHPRs). This term involves processes that may depend on the junctional SR [Ca2+] directly (as in the transition 
) or indirectly (as in the transition 
), as well as terms dependent on the membrane voltage (such as the transition 
). Using the decomposition of Q given by Eq. (17), one can see that [ρQ]i is a function of V and cjsr given by
where
provides the voltage-dependence, the superscripts of
, and
indicate row and column indices of these matrices,
ρj(
cjsr,
t) is the probability density for state
j, and
and
are given by Eqs.
(21).
The concentration balance equations governing the bulk myoplasmic (cmyo) and network SR (cnsr) [Ca2+] in the probability density formulation are identical to those used in the Monte Carlo approach (Eqs. (1)), except that the fluxes 
and 
are dependent on the densities 
that is,
where
is a function of
cjsr (Eq.
(21)).
Results
Moments of junctional SR [Ca2+]
The application of the moment-closure technique to the local control model of Ca2+-induced Ca2+ release (CICR) in cardiac myocytes presented above begins by writing the qth moment of the univariate probability density function, ρi(cjsr, t), as
where the nonnegative integer
q indicates the moment degree in
and is an exponent in (
cjsr)
q. As defined in Eq.
(18),
ρi(
cjsr,
t) is the distribution of [Ca
2+] in a large number of junctional SR compartments jointly distributed with CaRU state. Thus, the zero
th moment
corresponds to the probability—denoted as
πi(
t) in Williams et al.
18—that a randomly sampled CaRU is in state
i; that is,
where conservation of probability implies
Because the joint probability densities do not individually integrate to unity, the first moment,
is related to the expected value of the junctional SR [Ca
2+] conditioned on CaRU state through
while the conditional variance of the junctional SR [Ca
2+] is
Expressing fluxes in terms of moments
Considering Eqs. (1) and Eqs. (26), one sees that the fluxes 
and 
mediate the influence of the distribution of diadic subspace and junctional SR [Ca2+] on the dynamics of the bulk myoplasmic [Ca2+] (cmyo) and the network SR [Ca2+] (cnsr). Using the definition of the moments of junctional SR [Ca2+] (Eq. (28)), these fluxes become functions of the zeroth and first moments,
Similarly, the total flux through all the L-type Ca
2+ channels (
, Eq.
(24)) and the RyR Ca
2+ channels
become
and
Note that the average diadic subspace and junctional SR Ca
2+ concentrations can also be written in terms of the moments,
and
and
when expressed in terms using these quantities.
Derivation of moment equations
Differentiating Eq. (28) with respect to time and using the equations of the univariate probability density approach (Eqs. (19)), we obtain a system of ODEs that describe the time-evolution of these moments defined in Eq. (28),
where
M=12, 1≤
i≤
M,
q=0, 1, 2, …, and
and
are given by Eqs.
(22). In this expression the CaRU model is specified by the
M×
M matrices
Kϕ,
Kds, and
Kjsr defined in Eq.
(17), and the superscripts in
,
and
indicate the transition rate or bimolecular rate constant in the
jth row and
ith column of these matrices. Note that, in the
M equations for the zero
th moments
the first two terms evaluate to zero because
q=0. When
q≥1, the first term depends on both the network SR [Ca
2+] (
cnsr) and the bulk myoplasmic [Ca
2+] (
cmyo) through
. The terms in the first summation have a similar dependence on
cmyo and this can affect transitions mediated by diadic subspace Ca
2+
, and the magnitude of these terms depends also on voltage through
. Perhaps most importantly, the presence of diadic subspace and junction SR Ca
2+-mediated transitions in the CaRU model implies that
is a function of
whenever
or
is nonzero. That is, Eq.
(36) is an open system of the form
where we write
as a shorthand for
. Consequently, Eq.
(36) is unusable in its current form, because to determine the time-evolution of the
qth moments one needs to know the value of the (
q+1)
th moments.
Moment closure
To utilize Eq. (36), we truncate the open system at the second moment (q=2) and close the system of ODEs by assuming that the third moment can be expressed as an algebraic function ϕ of the lower moments 
; that is,
The remainder of this section derives the required expression of the form
(Eqs.
(48)). This is accomplished by specifying the function
ϕ in a manner that would be strictly correct if the probability density functions were scaled
β-distributions. Note that choosing this form of
ϕ to perform the moment closure given by Eqs.
(39) is not equivalent to assuming that the probability density functions are well approximated by
β-distributions. What we are assuming is that the relationship between
and the lower moments
is similar to the relationship observed in the
β-distribution. This assumption is validated a posteriori by evaluating the accuracy of results obtained using this approach (see
Figs. 2–6).
The derivation begins by considering a random variable 
that is functionally dependent on 
through
In this expression, the minimum and maximum values of junctional SR [Ca
2+] are given by
and
where
are the steady-state values of
cjsr found by setting
in Eq.
(20),
where
and
are given by Eqs.
(22). In this way, the maximum and minimum junctional SR Ca
2+ concentrations are determined to be
where
If the probability density for
conditioned on CaRU state
i were
β-distributed, then
where the
β-function
B(
αi,
βi) appears as a normalization constant and
,
αi(
t), and
βi(
t) are all functions of time. Under this assumption, the first several conditional moments of
would be
and inverting these expressions gives
Note that Eq.
(42) implies the following relationship between the conditional moments of
and
,
where we have used
for
q=0, 1, and 2; consequently,
αi and
βi can be found as a function of
and
. These parameters allow us to approximate the third conditional moment of
,
which, in turn, allows us to approximate the third conditional moment of junctional SR [Ca
2+] given by
,where
After some simplification one obtains
which is an expression that takes the form
as required by Eq.
(41), because
is a function of
, and
given by Eqs.
(48).
Note that the expression 
derived above is one of several possibilities that we tested, but the only one that could be validated. For example, when ϕ was chosen in a manner that would be strictly correct if the probability densities were scaled normal or log-normal distributions, the resulting moment closure did not perform well (not shown). Using the β-distribution to derive ϕ makes sense because it is a continuous distribution defined on a finite interval. In addition, for particular values of αi and βi, the β-distribution (while remaining integrable) diverges at the boundaries 
. Similarly, prior work has established that the densities ρi(cjsr, t) can accumulate probability at the minimum and maximum junctional SR Ca2+ concentrations (Eqs. (43)) and diverge as 
or 
18,22. As mentioned above, the use of the β-distribution to derive ϕ is ultimately validated by evaluating the accuracy of results obtained using this approach (see Figure 2 and Figure 3 and Figure 4 and Figure 5 and Figure 6).
Representative Monte Carlo and moment-closure results
Fig. 2 shows representative results from the minimal whole cell model of EC coupling described above. In this simulated voltage-clamp protocol, the holding potential of −80mV is followed by a 20-ms duration test potential to −10mV. The Monte Carlo result (shaded line) which involves a large but finite number of Ca2+ release units (N=1000) can be easily spotted by the fluctuations due to the stochastic gating of the CaRUs. The moment-closure result (solid line) that assumes N→∞ lacks these fluctuations. The top and bottom panels of Fig. 2 show the average diadic subspace 
and junctional SR 
Ca2+ concentrations in the Monte Carlo calculation (shaded lines) as well as the corresponding quantities from the moment-closure calculation (solid lines, Eqs. (34)). The middle two panels of Fig. 2 show the total Ca2+ influx through L-type Ca2+ channels 
and the total Ca2+ release from the RyR Ca2+ channels 
for the Monte Carlo calculation (shaded lines) as well as the corresponding quantities for the moment-closure result (solid lines, Eqs. (32)). In both the Monte Carlo and moment-closure calculations, the test potential of −10mV leads to 16× gain, here defined as the ratio 
where the overbar indicates an average over the duration of the pulse.
Fig. 3 shows [Ca2+] in the bulk myoplasm (cmyo) and network SR (cnsr) before, during, and after the −10mV voltage pulse (note change in timescale). In both cases the moment-closure result is shown as a solid line while the Monte Carlo is displayed as a dashed line (note agreement). While junctional SR depletion develops rapidly after the initiation of the voltage pulse (not shown), refilling the junctional SR compartments via diffusion of Ca2+ from the network SR (
in Eq. (2)) depletes this compartment (cnsr), which does not fully recover until ∼300ms after the termination of the voltage pulse.
Taken together, Figure 2 and Figure 3 validate our implementation of both the Monte Carlo and moment-closure approaches. Also note that the similarity of these results to Figure 2 and Figure 3 in Williams et al. 18 indicates that the six-state RyR megachannel model (Eq. (13))—used here because it takes the form of Eq. (17)has behavior similar to the two-state model of Williams et al. 18.
Dynamics of the moments of junctional SR [Ca2+]
The top row of Fig. 4 shows the time evolution of the probability of three selected CaRU states during the simulated voltage-clamp protocol of Figure 2 and Figure 3, as calculated using both the Monte Carlo (shaded lines) and moment-closure (solid lines) methods. Before the voltage pulse, the probability of state 
(DHPR in state 
and RyR in state 
, see Eqs. (13)) is ∼1, but during the voltage pulse to −10mV this probability drops to ∼0.78 (20–40ms). Conversely, the probability of CaRU state 
(DHPR open and RyR open) and 
(DHPR closed and RyR open) both increase during the voltage pulse. The dynamics of voltage-dependent activation of DHPRs and subsequent triggering of the opening of RyR megachannels is similar in both the Monte Carlo (shaded lines) and moment-closure (solid lines) calculations.
The second row of Fig. 4 shows the mean junctional SR [Ca2+] conditioned on CaRU state for the Monte Carlo (shaded line) and the moment-closure (solid line) techniques. In the Monte Carlo calculation this conditional mean is given by
where
Ni(
t) is the number of CaRUs in state
i at time
t and
so that the sum includes only those CaRUs in state
i. The corresponding quantity in the moment-closure technique is the conditional expectation
(Eq.
(29)). Note that before the voltage pulse the expectation of SR [Ca
2+] is ∼1000
μM when conditioning on CaRU state
, 851
μM when conditioning on CaRU state
, and 306
μM when conditioning on CaRU state
. That is, at the holding potential of −80mV, the stochastic gating of CaRUs leads to depletion of junctional SR [Ca
2+] associated with release sites with open RyR megachannels (more pronounced in
than
because the former state is longer-lived). However, the probability of CaRU states
and
is very low at −80mV and, consequently, the expectation of junctional SR [Ca
2+] irrespective of CaRU state given by the weighted average
in the Monte Carlo model and
in the moment-closure calculation is ∼1000
μM, consistent with
Fig. 2. Also note that during the voltage pulse the conditional expectation of junctional SR [Ca
2+] decreases for CaRU states
and
but first increases and then decreases for CaRU state
presumably because the increasing probability of state
during the pulse is due to CaRU transitions into this state from others (such as
) that have higher resting junctional SR [Ca
2+].
The third row of Fig. 4 shows the variance of the junctional SR [Ca2+] conditioned upon the CaRU state for the Monte Carlo (shaded line) and the moment-closure (solid line) techniques. For the Monte Carlo calculation
where
Ni and
ni(
t) are defined as in Eq.
(54), while the corresponding conditional variance of the junctional SR [Ca
2+] in the moment-closure calculation is
(Eq.
(30)). Note that during the voltage pulse the conditional variance of
cjsr increases, as the dynamics of EC coupling lead to increased heterogeneity of junctional SR [Ca
2+], and that the moment-closure technique accurately accounts for this heterogeneity (compare
shaded and
solid lines).
The distribution of junctional SR [Ca2+] conditioned on CaRU state
Fig. 5 shows a snapshot of the distribution of junctional SR [Ca2+] (cjsr) conditioned upon the state of the Ca2+ release unit at t=30ms, midway through the voltage pulse protocol of Figure 2 and Figure 3 and Figure 4. For clarity, the five closed states of the RyR megachannel (
in Eq. (13)) have been lumped resulting in a contracted presentation with four CaRU states: 
, and 
where 
and 
(Eq. (16)). Thus, the two histograms on the bottom of Fig. 5 indicate the distribution of JSR [Ca2+] when the DHPR is open 
, while the two histograms on the right of Fig. 5 indicate the distribution of JSR [Ca2+] when the RyR megachannel is open 
.
Fig. 5 shows a broad range of junctional SR [Ca2+] regardless of CaRU state, consistent with the high variances at t=30ms in Fig. 4. For example, when the RyR megachannel is closed (
and 
, left panels), a randomly sampled junctional SR is likely to be replete, as indicated by the large vertical bar at cjsr≈1000μM. However, one can also find depleted junctional SR [Ca2+] associated with closed RyR megachannels, where RyRs have recently opened and the junctional SR has not had time to refill. Conversely, when the RyR is open (
and 
right panels), the probability mass has shifted to lower junctional SR [Ca2+].
The diamonds of Fig. 5 show β-distributions with the same mean and variance as the histograms obtained from Monte Carlo simulation. While the agreement is noteworthy, this correspondence is not required for the moment-closure technique to work well. What is required is that the relationship between the third 
and lower 
moments in the histograms is similar to that observed in the β-distribution. For example, the histogram junctional SR [Ca2+] for CaRU state 
at t=30ms has moments of 
μM, 
, and 
. When moments 0–2 are used to estimate the third moment using Eq. (53) with 
and 
(Eqs. (43)), one obtains 
, for a relative error of only 0.1%. It is this low error that is responsible for the excellent agreement between the moment-closure result and the Monte Carlo calculation observed in Figure 2 and Figure 3 and Figure 4.
The model displays gain and gradedness
To further validate the moment-closure approach by comparison to Monte Carlo simulation, Fig. 6 summarizes a large number of simulated whole-cell voltage-clamp protocols such as those presented in Figure 2 and Figure 3 and Figure 4. The open circles of Fig. 6 show the trigger Ca2+ influx via L-type Ca2+ channels integrated over the 20-ms voltage step to test potentials in the range −40 to 40mV using 1000 CaRUs (the plot is normalized to maximum value of 
). The dotted line of Fig. 6 shows that the trigger Ca2+ influx in the moment-closure calculation agrees with the Monte Carlo simulations. Similarly, the open squares of Fig. 6 show the voltage-dependence of the Ca2+ release flux (normalized to maximum value of 
), while the dashed lines of Fig. 6 show that the Ca2+ release flux observed in the moment-closure calculation agrees with the Monte Carlo simulations. Note that the Monte Carlo and moment-closure calculations exhibit graded Ca2+ release. Furthermore, the EC coupling gain 
is a decreasing function of voltage, in the range of 32–15× for test potentials between −40 and 0mV. Most importantly, the Monte Carlo and moment-closure calculations are nearly identical (compare open diamonds and solid line).
Computational efficiency of the moment-closure approach
While the previous sections have shown that the moment-closure and Monte Carlo calculations are essentially equivalent in terms of the dynamic cellular responses they predict, it is important to note that the moment-closure approach is significantly faster than Monte Carlo simulation. The Monte Carlo simulations presented above are performed using Δt=0.01 μs, a value chosen so the probability of transition occurring in each CaRU is <5% per time step. Table 1 shows that the run time for these 60-ms simulations increases approximately linearly with the number of CaRU units; for example, an N=10,000 simulation takes ∼11 times longer than a N=1000 simulation. When our current implementation of the moment-closure method is employed using a nonadaptive time step of Δt=0.01 μs, the run time is 95min, which is ∼100 times faster than Monte Carlo simulations with a physiologically realistic number of CaRUs (e.g., N=10,000). However, a time step of 0.01 μs is much smaller than required for integrating the moment-closure ODEs. When this artificial constraint is removed and the moment-closure approach is benchmarked using a nonadaptive time step as large as numerical stability will allow, the calculations are 8755:0.9=9728 times faster than Monte Carlo simulations containing N=10,000 CaRUs. That is, the computational efficiency of the moment-closure approach is nearly four orders-of-magnitude superior to physiologically realistic Monte Carlo simulations, while leading to nearly identical results (see Figure 2 and Figure 3 and Figure 4 and Figure 6). Furthermore, integration methods that utilize adaptive time-stepping are likely to further enhance the computational advantage of the moment-closure approach to modeling local control of EC coupling.
, Jin, Jncx, Jserca, and Jleak are described in the text and the Appendix .
, total Ca2+ flux via the DHPR Ca2+ channels 
, total Ca2+-induced Ca2+ release flux 
, and average junctional SR [Ca2+] 
. The Monte Carlo simulation used N=1000 Ca2+ release units and parameters as in Table 2,Table 3,Table 4.
, and the conditional variance of cjsr
, for three selected CaRU states, 
(left column), 
(middle column), and 
(right column). The Monte Carlo simulation used N=2000 Ca2+ release units.
, DHPR and RyR megachannel both closed; 
, DHPR open and RyR megachannel closed; 
, DHPR closed and RyR megachannel open; and 
, DHPR and RyR megachannel both open.
is shown as open circles (Monte Carlo) and dotted line (moment closure). Integrated RyR flux 
is shown as open squares (Monte Carlo) and dashed line (moment-closure). EC coupling gain (
, right axis) is shown as open diamonds (Monte Carlo) and solid line (moment-closure).