## Abstract

Population density approaches to modeling local control of Ca^{2+}-induced Ca^{2+} release in cardiac myocytes can be used to construct minimal whole cell models that accurately represent heterogeneous local Ca^{2+} signals. Unfortunately, the computational complexity of such “local/global” whole cell models scales with the number of Ca^{2+} release unit (CaRU) states, which is a rapidly increasing function of the number of ryanodine receptors (RyRs) per CaRU. Here we present an alternative approach based on a Langevin description of the collective gating of RyRs coupled by local Ca^{2+} concentration ([Ca^{2+}]). The computational efficiency of this approach no longer depends on the number of RyRs per CaRU. When the RyR model is minimal, Langevin equations may be replaced by a single Fokker-Planck equation, yielding an extremely compact and efficient local/global whole cell model that reproduces and helps interpret recent experiments that investigate Ca^{2+} homeostasis in permeabilized ventricular myocytes. Our calculations show that elevated myoplasmic [Ca^{2+}] promotes elevated network sarcoplasmic reticulum (SR) [Ca^{2+}] via SR Ca^{2+}-ATPase-mediated Ca^{2+} uptake. However, elevated myoplasmic [Ca^{2+}] may also activate RyRs and promote stochastic SR Ca^{2+} release, which can in turn decrease SR [Ca^{2+}]. Increasing myoplasmic [Ca^{2+}] results in an exponential increase in spark-mediated release and a linear increase in nonspark-mediated release, consistent with recent experiments. The model exhibits two steady-state release fluxes for the same network SR [Ca^{2+}] depending on whether myoplasmic [Ca^{2+}] is low or high. In the later case, spontaneous release decreases SR [Ca^{2+}] in a manner that maintains robust Ca^{2+} sparks.

- Langevin equation
- Fokker-Planck equation
- calcium release site
- multiscale whole cell model
- calcium homeostasis

intracellular calcium (Ca^{2+}) signaling involves a complex interplay between global (cell-wide) changes in Ca^{2+} concentration ([Ca^{2+}]) and local (subcellular) Ca^{2+} release events. Local signals are caused by plasma membrane Ca^{2+} influx and release of Ca^{2+} from intracellular stores, primarily the endoplasmic/sarcoplasmic reticulum (ER/SR). Spatially localized Ca^{2+} release events mediated by clusters of intracellular Ca^{2+} channels, IP_{3} receptors (IP_{3}Rs) or ryanodine receptors (RyRs) on the ER/SR membrane, are referred to as “Ca^{2+} sparks” or “puffs” (see Ref. 2 for review).

While plasma membrane ion channels in a small cell experience essentially the same time course of membrane voltage, intracellular Ca^{2+} channels experience radically different local [Ca^{2+}], even during global Ca^{2+} responses, and clusters of IP_{3}Rs and RyRs are in fact only locally coupled via the buffered diffusion of intracellular Ca^{2+}. That is, when one or several of the channels in a Ca^{2+} release unit (CaRU) are open, the [Ca^{2+}] experienced by spatially localized channels is dramatically different from the [Ca^{2+}] in the bulk myoplasm. For this reason, conventional whole cell modeling of Ca^{2+} dynamics based on Hodgkin-Huxley-like gating variables for the dynamics of intracellular channels is not always appropriate.

Mechanistic models of ER/SR Ca^{2+} release often represent the stochastic gating of Ca^{2+} channels using Monte Carlo methods. When these approaches are applied to cardiac myocytes, voltage-gated L-type Ca^{2+} channel(s) interact with a cluster of RyRs through changes in [Ca^{2+}] in small “dyadic subspaces” between the sarcolemmal and SR membranes. These models also sometimes consider depletion of junctional SR [Ca^{2+}] that may influence Ca^{2+} spark termination and refractoriness (31, 32, 35). Realistic global (cell-wide) SR Ca^{2+} release can be reproduced by Monte Carlo simulation of the stochastic triggering of sparks from hundreds to thousands of CaRUs (19, 20, 29, 32). However, such simulations of local control of excitation-contraction coupling are computationally demanding, especially when each CaRU is composed of interacting Markov chain models of individual RyRs (e.g., see Ref. 23).

Population density approaches are an alternative to Monte Carlo simulations that produce realistic and computationally efficient models by using a master equation to represent heterogeneous local Ca^{2+} signals in dyadic subspaces and junctional SR domains (37). This approach involves the numerical solution of advection-reaction equations for the time-dependent bivariate probability density of subspace and junctional SR [Ca^{2+}] conditioned on CaRU state, coupled to ordinary differential equations (ODEs) for the bulk myoplasmic and network SR [Ca^{2+}]. This methodology was validated in prior work (37) and an associated moment-based approach to simulating the probability distribution of junctional SR [Ca^{2+}] was benchmarked to be several orders of magnitude faster than conventional Monte Carlo simulation of the dynamics of local Ca^{2+} associated with a physiological number of CaRUs (38).

One disadvantage of the population density approaches to modeling local control is that their run times (computational efficiency) are proportional to the number of CaRU states. When realistically modeled as the collective gating of identical and indistinguishable RyRs, the number of CaRU states is exponential in the number of channel states. Population density and moment-based methods for multiscale (i.e., local/global) whole cell modeling are limited by this state-space explosion.

Here we present an alternative local/global whole cell modeling approach based on a Langevin formulation of the stochastic Ca^{2+} release via CaRUs. We assume that the number of RyRs per CaRU is large enough that the fraction of channels in each state can be treated as a continuous variable. We show that the Langevin description of the collective gating of RyRs is a good approximation to the corresponding discrete-state continuous-time Markov chain model when the number of RyRs per release site is in the physiological range. By coupling the numerical solution of such Langevin equations to balance equations for the bulk myoplasmic and network SR [Ca^{2+}], a local/global whole cell model is produced whose run time scales with the number of states in the Markov chain model for an individual RyR, as opposed to the far greater number of states in a compositionally defined CaRU. When the RyR model is minimal, these Langevin equations may be replaced by a single Fokker-Planck equation for a randomly sampled CaRU, yielding an extremely compact and efficient local/global whole cell model. We illustrate the usefulness and computational efficiency of the Fokker-Planck equation-based local/global whole cell model by performing parameter studies motivated by recent experiments (3, 40).

In intact ventricular myocytes of the healthy heart, the balance of diastolic SR Ca^{2+} leak and uptake maintains the appropriate SR Ca^{2+} load. While the SR Ca^{2+} leak is mediated primarily by RyRs, the contributions of spark- and nonspark-mediated SR Ca^{2+} release depend on the concentration of both myoplasmic and SR [Ca^{2+}] (10, 27, 30, 40). When SR [Ca^{2+}] is low, SR Ca^{2+} leak occurs primarily through spark-independent pathways. Conversely, when SR [Ca^{2+}] is high, spontaneous Ca^{2+} sparks make a large contribution to SR leak. In pathophysiological conditions that include SR Ca^{2+} overload, increased SR Ca^{2+} leak may generate spontaneous sparks that trigger Ca^{2+}-induced Ca^{2+} release (CICR) from neighboring CaRUs, thereby initiating arrhythmogenic spontaneous Ca^{2+} waves (9).

With the use of permeabilized ventricular myocytes, a reduced experimental preparation that allows precise control of myoplasic [Ca^{2+}], Bovo et al. (3) observed that increasing myoplasmic [Ca^{2+}] results in an exponential increase in spark-mediated release and a linear increase in nonspark-mediated release. These results are reproduced by the Fokker-Planck equation-based local/global whole cell model that is the focus of this article. In addition, the model predicts potentially significant characteristics of Ca^{2+} homeostasis in permeabilized cells. For example, in the local/global whole cell model, two distinct steady states may exist for a given network SR [Ca^{2+}]. One steady-state corresponds to low myoplasmic [Ca^{2+}] and small SR Ca^{2+} release flux that is dominated by stochastic leak, while the other corresponds to high myoplasmic [Ca^{2+}] and large release flux mediated by Ca^{2+} sparks. Interestingly, for any clamped myoplasmic [Ca^{2+}] that is large enough to trigger spark-mediated release, the local/global model predicts that the resulting spontaneous stochastic Ca^{2+} release tends to decrease the network SR Ca^{2+} load just enough to maintain robust Ca^{2+} sparks.

## METHODS

#### Markov chain description of a Ca^{2+} release site.

The most straightforward starting point for the presentation of the Langevin description of a Ca^{2+} release site (CaRU) is the following two-state Markov chain model of a stochastically gating RyR,
(1)

where c is the local [Ca^{2+}], *k*^{+}c^{η} and *k*^{−} are transition rates with units of reciprocal time, *k*^{+} is an association rate constant with units of concentration^{−η}·time^{−1}, and η is the cooperativity of Ca^{2+} binding. Under the assumption that a collection of *N* two-state RyRs are instantaneously coupled by a local [Ca^{2}] associated with the RyR cluster, the transition diagram for the CaRU as a collective entity is (8)
(2)

where the states {0, 1, . . ., *N*} correspond to the number of open channels (*N*_{O}) and c_{n} is the local [Ca^{2+}] experienced by RyRs when *N*_{O} = *n*.

Figure 1*A* shows a Markov chain simulation of a CaRU composed of *N* = 20 two-state channels. For simplicity we here assume that the local [Ca^{2+}] is a linear function of *N*_{O}, that is,
(3)

where c_{∞} is the bulk or background [Ca^{2+}] and c_{*} determines the increment in local [Ca^{2+}] following an individual RyR opening. The corresponding relationship between *N*_{O} and local [Ca^{2+}] is more realistic in the local/global whole cell model (*Eqs. 24* and *25*).

#### Langevin Ca^{2+} release site model.

We will write *f*_{O}(*t*) as the time-varying fraction of open RyRs, that is,
(4)

The Langevin equation that corresponds to a CaRU composed of *N* two-state channels (see above) is a stochastic ordinary differential equation (SDE) of the form
(5)

where c̄ = *N*c_{*} and ξ(*t*) is a rapidly varying forcing term (Gaussian white noise) with zero mean
(6)

As discussed in appendix a, the magnitude of the noise term, ξ(*t*), is characterized by the two-time covariance (16, 26),
(7)

where δ is the Dirac delta function and γ(*f*_{O}) is the infinitesimal variance of *f*_{O} and is given by
(8)

With the use of parameters that lead to Ca^{2+} sparks, Fig. 1*B* shows that the Langevin simulation of a 20-channel CaRU is qualitatively similar to the corresponding Markov chain simulation (Fig. 1*A*).

#### Equivalence of Markov chain and Langevin formulations.

The Langevin CaRU model is expected to well approximate the Markov chain model when the number of RyRs per CaRU (*N*) is sufficiently large. To determine whether this convergence occurs for a physiological number of RyRs (10–200 per CaRU in skeletal and cardiac myocytes; Ref. 11), we compare the stationary distributions for *N*_{O} (Fig. 1*C*). The white histogram of Fig. 1*C* shows that the bimodal Markov chain stationary distribution has a local maxima near *N*_{O} = 0 and *N*_{O} = 15. This distribution reflects the fact that the *N* RyRs are usually closed but occasionally open in concert as is characteristic for Ca^{2+} sparks. Comparison to the corresponding distribution of the Langevin model (black histogram) shows that the SDE formulation is a good approximation to the Markov chain, even when the number of RyRs per CaRU is on the low end of the physiological range.

In the Langevin CaRU formulation, the state space for *f*_{O} is continuous (0 ≤ *f*_{O} ≤ 1). The Fokker-Planck equation solved by the probability density function for the fraction of open channels, ρ(*f*, *t*), is given by (12)
(9)

where ρ(*f*,*t*)d*f* = Pr{*f* ≤ *f*_{O}(*t*) < *f* + d*f*}. Note that in these expressions, *f*_{O} is the random variable and *f* is the independent variable of the probability density. The drift and diffusion terms in *Eq. 9* are given by
(10)
(11)

where *v*^{±}(*f*) are the rates of the elementary processes leading to an increase and decrease in the fraction of open channels, that is,
(12)
(13)

and c̄ = c_{*}*N* as above (*Eq. 5*).

Setting the left-hand side of *Eq. 9* equal to zero ( = 0), denoting the stationary density by ρ_{ss}(*f*) and applying boundary conditions ρ_{ss}(*f*) → 0 as *f* → ±∞, it can be shown that (12)
(14)

where θ is a normalization constant such that ∫ρ_{ss}(*f*)d*f* = 1 and
(15)

is an accumulation function with a lower limit of integration satisfying α(*f*)/γ(*f*) = 0 for *f* ≤ *a*. In fact, *U* may be any antiderivative satisfying *U*′ = α/γ, because the normalization of ρ_{ss} determines the constant of integration.

Figure 1*D* shows the stationary density ρ_{ss}(*f*) for the 20-channel Fokker-Planck CaRU model described above. The + symbols in Fig. 1*C* are binned values of ρ_{ss}(*f*) that may be compared with (and agree with) the stationary distributions of the Markov chain (white histogram) and Langevin (black histogram) descriptions. appendix b provides more comparisons of Markov chain, Langevin, and Fokker-Planck CaRU simulations.

#### Full local/global whole cell model.

Having validated the Langevin CaRU model in the previous sections, we are prepared to construct the local/global whole cell model of Ca^{2+} homeostasis in permeabilized ventricular myocytes that is the focus of this article. Figure 2 shows the relationship between the bulk Ca^{2+} concentrations of the myoplasm (c_{myo}) and the network SR (c_{nsr}) and the local Ca^{2+} concentrations associated with each CaRU. With respect to global aspects of Ca^{2+} signaling, the material balance equations of the whole cell model are
(16)
(17)

where Λ_{nsr} is an effective volume ratio that accounts for both physical volume and Ca^{2+} buffering capacity of the myoplasm and network SR. A plasma membrane flux may take the form *J*_{pm} = *k*_{pm}(c_{ext} − c_{myo}). The sarco(endo)plasmic reticulum Ca^{2+}-ATPase (SERCA) type Ca^{2+} ATPase flux is (37)
(18)

The aggregate fluxes *J*_{myo}^{T} = ∑_{m = 1}^{M}*J*_{myo}^{m} and *J*_{nsr}^{T} = ∑_{m = 1}^{M}*J*_{nsr}^{m} in *Eqs. 16* and *17* account for the stochastic dynamics of Ca^{2+} release, where *J*_{myo}^{T} = *v*_{myo}(c_{ds}^{m} − c_{myo}) with *v*_{myo} = *v*_{myo}^{T}/*M* is the flux from the *m*th dyadic subspace into the bulk myoplasm and *J*_{nsr}^{m} = *v*_{nsr}(c_{nsr} − c_{jsr}^{m}) where *v*_{nsr} = *v*_{nsr}^{T}/*M* is the flux from the network SR to the *m*th junctional SR (*m* = 1, 2, . . ., *M*). See Table 1 for parameters.

Each CaRU in the whole cell model is a collection of *N* RyRs with open fraction *f*_{O}^{m} and associated dyadic subspace (c_{ds}^{m}) and junctional SR (c_{jsr}^{m}) Ca^{2+} concentrations:
(19)
(20)
(21)

In these equations, λ_{ds} and λ_{jsr} are effective volume ratios, that is, λ_{ds} = (*V*_{ds}/β_{ds})/(*V*_{myo}/β_{myo}), where *V*_{ds} = *V*_{ds}^{T}/*M* and *V*_{ds}^{T} is the aggregate volume of the diadic subspaces (similarly for λ_{jsr}). *J*_{rel}^{m} is the release flux though the *m*th RyR cluster given by *J*_{rel}^{m} = *v*_{rel}*f*_{O}^{m}(c_{jsr}^{m} − c_{ds}^{m}) for *m* = 1, 2, . . ., *M* and *v*_{rel}^{T}/*M*. The random functions of time ξ^{m}(*t*) are indepenent Gaussian white noise terms with zero mean, = 0 for all *m*, and the two-time covariances are
(22)

where (23)

Note that the dyadic subspaces only influence each other through the bulk concentrations c_{myo} and c_{nsr}. Below we refer to *Eqs. 16–23* as the “full local/global whole cell model.”

#### Reduced local/global whole cell model.

The Langevin description of each CaRU (*Eqs. 19–21*) in the full local/global model may be simplified by assuming that the dyadic subspace and junctional SR rapidly equilibrate with the bulk myoplasmic and network SR [Ca^{2+}], that is, *J*_{rel}^{m} = *J*_{myo}^{m} and *J*_{rel}^{m} = *J*_{nsr}^{m}. These balanced fluxes relate the 2*M* domain Ca^{2+} concentrations, c_{ds}^{m} and c_{jsr}^{m}, to the bulk concentrations, c_{myo} and c_{nsr}, and the fraction of open channels *f*_{O}^{m} in the *m*th CaRU as follows (22),
(24)
(25)

where χ_{myo}^{m} = /(*v*_{myo} + ) and χ_{nsr}^{m} = /(*v*_{nsr} + ), *v*_{nsr}/( + *v*_{nsr}), *v*_{myo}/( + *v*_{myo}), = *v*_{rel} *f*_{O}^{m}. *Equations 24* and *25* eliminate 2*M* of the 3*M* ODEs representing the population of *M* CaRUs, with the remaining ODEs,
(26)

dependent on the rapidly equilibrated dyadic subspace concentration c̄_{m}^{ds} that is an algebraic function of *f*_{O}^{m}, c_{myo}, and c_{nsr}. Realizations of this “reduced local/global whole cell model” are obtained by numerically integrating *Eqs. 16*, *17*, and *26*.

#### Fokker-Planck local/global whole cell model.

The full and reduced local/global whole cell models presented above include heterogeneous local Ca^{2+} signaling and stochastic Ca^{2+} release. Unfortunately, a physiologically realistic ventricular myocyte simulation would involve *M* ≈ 20,000 CaRUs (6). Rather than perform Monte Carlo simulations with a lesser, unphysiological value for *M* that is computationally feasible, we recognize that a Fokker-Planck equation similar to *Eq. 9* is the master equation for a CaRU and its associated domains. Because the *M* CaRUs in the whole cell model are identical and independent except for fluxes to and from the bulk myoplasm and network SR, we replace the *M* SDEs representing these CaRUs (*Eq. 26*) with this Fokker-Planck equation (a good approximation for large *M* that is exact as *M* → ∞). In this way, we obtain the “Fokker-Planck local/global whole cell model.”

In the study of Ca^{2+} homeostasis in permeabilized ventricular myocytes presented below, the governing equations are *Eqs. 9*, *16*, and *17*, with the fluxes *J*_{myo}^{T} and *J*_{nsr}^{T} redefined as functions of the probability distribution function for *f*_{O} in a randomly sampled CaRU. In permeabilized myocytes, the bulk myoplasmic [Ca^{2+}] is clamped (*k*_{pm} is large) and c_{myo} ≈ c_{ext} is no longer a variable but a parameter. Consequently, *Eq. 16* is superfluous and the governing equations for the Fokker-Planck equation description of the local/global model of permeabilized ventricular myocytes are therefore given by
(27)
(28)

where (29)

In *Eq. 28*, α(*f*) and γ(*f*) are given by *Eqs. 10* and *11*, with
(30)
(31)

The equilibrated domain concentrations are given by (32) (33)

where χ_{myo} and χ_{nsr} are the following functions of *f*,
(34)
(35)

where (36) (37)

and = *v*_{rel}^{T} (cf. *Eqs. 24* and *25*). In the local/global whole cell model calculations presented below, the Fokker-Planck equation was numerically integrated using a total variation diminishing scheme (37), with boundary conditions as described in appendix c.

## RESULTS

#### Calcium homeostasis in the local/global whole cell model.

We use the Fokker-Planck version of the reduced local/global model (*Eqs. 27–33*) to investigate Ca^{2+} homeostasis in permeabilized ventricular myocytes, in particular, the influence of c_{myo} on SR Ca^{2+} load and release. The relationship between c_{myo} and Ca^{2+} homeostasis is complex, as c_{myo} can promote elevated c_{nsr} through increased SERCA uptake. On the other hand, a sufficiently elevated c_{nsr} also promotes Ca^{2+} sparks that may deplete the network SR (i.e., decrease c_{nsr}).

With the use of an intermediate value for the myoplasmic [Ca^{2+}] (c_{myo} = 0.18 μM) in the permeabilized ventricular myocyte model, Fig. 3*A* shows the bimodal steady-state probability density function for the fraction of open channels, ρ_{ss}(*f*), calculated via the Fokker-Planck version of the whole cell model (solid line). This bimodal density reflects the dynamics of CaRUs composed of RyRs that are usually closed but occasionally open in a concerted fashion. For comparison, Fig. 3*A* also shows a (nearly identical) estimate of the steady-state density function obtained from a whole cell model with the corresponding Langevin description of *M* = 200 release sites (dashed curve). Figure 3*B* compares the stationary distribution for a whole cell model that uses a Markov chain description of release sites (white bars) and the corresponding distribution calculated via the Fokker-Planck version of the model (appropriately discretized). The two histograms are qualitatively identical and in strong qualitative agreement (the Markov chain simulation is slightly shifted to larger *f*_{O}), validating the minimal Fokker-Planck formulation of the local/global whole cell model (*Eqs. 27–37*).

With the use of the Fokker-Planck-based whole cell model, Fig. 3*C* shows the monotone increasing relationship between the fraction of open channels and stochastic Ca^{2+} release rate, given by *v*_{rel}*f*(c̄_{jsr} − c̄_{ds}) where c̄_{jsr} and c̄_{ds} are functions of *f* (*Eqs. 32* and *33*). Figure 3*D* shows the steady-state release flux density, given by *v*_{rel}*f*(c̄_{jsr} − c̄_{ds})ρ_{ss}, that is, the product of the curves in Fig. 3, *A* and *C*. Note that the steady-state release flux density is also a bimodal function of *f*, with the first and second modes corresponding to nonspark-mediated (light gray area, *f* < 0.1) and spark-mediated stochastic Ca^{2+} release (dark gray area, *f* ≥ 0.1), respectively.

Figure 4 shows steady-state values for total release flux (*J*_{rel}^{T}), network SR [Ca^{2+}] (c_{nsr}), and the spark *Score* as a function of c_{myo}, obtained from simulation of the local/global model using the Langevin (+ symbols) and Fokker-Planck (solid lines) descriptions of the CaRU population. The spark *Score* is the index of dispersion of the fraction of open channels (*f*_{O}),
(38)

where E[*f*_{O}] = ∫*f*ρ_{ss}d*f*, Var[*f*_{O}] = ∫(*f* − E[*f*_{O}])^{2}ρ_{ss}d*f*, and ρ_{ss}(*f*) is the steady-state probability density of open channels. The spark *Score* takes values between 0 and 1, and a *Score* greater than ∼0.25 indicates the presence of robust Ca^{2+} sparks (21). Over a wide range of c_{myo} values, there is agreement among *J*_{rel}^{T}, c_{nsr}, and the spark *Score* calculated using the Langevin and Fokker-Planck approaches, validating the use of the Fokker-Planck version of the model and our implementation of both methods. Note that *J*_{rel}^{T} is a monotone increasing function of c_{myo} (Fig. 4*A*), while c_{nsr} is biphasic, increasing for c_{myo} < 0.2 μM and decreasing for c_{myo} > 0.2 μM (Fig. 4*B*). The spark *Score* shows similar biphasic dependence on c_{myo} (Fig. 4*C*).

The biphasic dependence of c_{nsr} and the spark *Score* on c_{myo} can be understood by considering the representative stochastic trajectories for the fraction of open channels in a randomly sampled CaRU in the Langevin model (Fig. 4*A*) or, alternatively, the steady-state population density function ρ_{ss}(*f*) in the Fokker-Planck model (Fig. 4*C*). For a low myoplasmic [Ca^{2+}] (c_{myo} = 0.1 μM), ρ_{ss}(*f*) is located near *f* = 0, consistent with few channel openings (Fig. 4, *A* and *C*, *insets*). As c_{myo} increases to an intermediate value of 0.2 μM, increased SERCA uptake elevates c_{nsr} and ρ_{ss}(*f*) is distinctly bimodal, consistent with robust sparks and the observed increase in *J*_{rel}^{T} and *Score*. However, a further increase in myoplasmic [Ca^{2+}] (c_{myo} = 0.6 μM) promotes tonic activation of CaRUs (as opposed to sparks, for a decreasing *Score*). The resulting increase in release flux (*J*_{rel}^{T}) depletes the network SR [Ca^{2+}] (lower values of c_{nsr}) and eliminates robust sparks. appendix d provides more details regarding the influences of c_{myo} on steady-state spark statistics.

In Fig. 5, *J*_{rel}^{T} (black dashed lines) and *J*_{pump} (solid gray lines) are shown as a function of c_{nsr} for three values of c_{myo}. *J*_{rel}^{T} is a monotone increasing function of c_{nsr}; the increasing slope at high c_{nsr} levels is due to spark-mediated Ca^{2+} release. *J*_{pump} decreases approximately linearly with c_{nsr}, and both *J*_{rel}^{T} and *J*_{pump} increase for increasing c_{myo}. The intersection of the *J*_{rel}^{T} and *J*_{pump} curves (open circles) indicate the steady-state total release flux and SR Ca^{2+} load (c_{nsr}) for a given value of c_{myo} (solid black line, arrow indicates increasing c_{myo}). For a given c_{nsr}, two distinct steady-states are possible, one with low c_{myo} and *J*_{rel}^{T} (primarily nonspark-mediated release) and another with high c_{myo} and *J*_{rel}^{T} (primarily spark-mediated release). The next section further explores the dependence of spark- and nonspark-mediated release on c_{myo}.

#### Spark- and nonspark-mediated SR Ca^{2+} release.

In a recent experimental study, Bovo et al. (3) demonstrated that myoplasmic Ca^{2+} levels augment both spark-mediated SR Ca^{2+} release and nonspark-mediated SR Ca^{2+} release in ventricular myocytes (3). While controlling myoplasmic [Ca^{2+}] (c_{myo}) by permeabilization of the cell plasma membrane, the time course of network SR [Ca^{2+}] (c_{nsr}) depletion was measured following application of the SERCA inhibitor thapsigargin (cf. Ref. 3, Fig. 1*A*). Assuming negligible SERCA activity (i.e., *J*_{pump} = 0), the rate of change of c_{nsr} was used as a measure of the SR Ca^{2+} release flux (see *Eq. 27*), and further analysis was performed to distinguish spark- and nonspark-mediated release as functions of c_{myo} and c_{nsr}. Figure 6 uses a similar protocol (setting *v*_{pump} = 0) to elucidate the influence of c_{myo} on spark- and nonspark-mediated release. Consistent with Bovo et al. and Fig. 4*B*, Fig. 6 shows that steady-state c_{nsr} increases as c_{myo} increases from 0.12 to 0.18 μM (compare initial values, solid, dashed, and thick solid lines). Consistent with the experiment, increasing c_{myo} in this range of concentrations also leads to an increased Ca^{2+} release rate, as evidenced by faster SR depletion upon simulated block of SERCA with thapsigargin (TG in Fig. 6).

Figure 7 shows the total release flux, *J*_{rel}^{T}, and the spark- and nonspark-mediated release (*J*_{rel}^{S} and *J*_{rel}^{NS} as defined in Fig. 3*C*) as a function of c_{nsr} during the SR depletion simulation of Fig. 6 for different values of c_{myo} (cf. Ref. 3, Fig. 3). While *J*_{rel}^{T}increases as a function of both c_{nsr} and c_{myo} (Fig. 7*A*), the contributions of the spark- and nonspark-mediated release (*J*_{rel}^{S} and *J*_{rel}^{NS}) are highly dependent on c_{nsr}. At low network SR [Ca^{2+}] (c_{nsr}), the spark-mediated release flux (*J*_{rel}^{S}) is negligible, but it increases exponentially as c_{nsr} increases (Fig. 7*B*). The nonspark-mediate release (*J*_{rel}^{NS}) is small for low c_{nsr} levels and increases as a linear function of c_{nsr} (Fig. 7*C*). When the SR load is clamped at c_{nsr} = 950 μM, both *J*_{rel}^{S} and *J*_{rel}^{NS} increase as c_{myo} increases (Fig. 7*D*). However, the spark-mediated release flux (*J*_{rel}^{S}) increases to a greater extent than the nonspark-mediated release (*J*_{rel}^{NS}). Steady-state calculations of *J*_{rel}^{T}, *J*_{rel}^{S}, and *J*_{rel}^{NS} closely agree with time-varying simulations (Fig. 7, + symbols). In summary, when the SR is depleted, most SR Ca^{2+} release occurs via nonspark-mediated release; conversely, when the SR is replete, most SR Ca^{2+} release occurs via Ca^{2+} sparks, more so as c_{myo} increases.

The number of RyRs per CaRU, *N*, can vary over a wide physiological range (11). Figure 8 shows the steady-state values for *J*_{rel}^{T}, *J*_{rel}^{S}, and *J*_{rel}^{NS} for different values of *N*. As *N* increases (scaling *v*_{rel} appropriately such that *J*_{rel}^{T} when all *N* channel are open is unchanged), *J*_{rel}^{T} becomes a steeper function of c_{nsr} (Fig. 8*A*). Interestingly, when the network SR [Ca^{2+}] is higher (c_{nsr} = 1,000 µM), *J*_{rel}^{T} is larger for large *N*, but when c_{nsr} is slightly lower (c_{nsr} = 950 μM), *J*_{rel}^{T} is smaller for large *N* (Fig. 8*A*, arrows). Spark-mediated release (*J*_{rel}^{S}) varies with *N* in a manner similar to *J*_{rel}^{T} (Fig. 8*B*), while nonspark-mediated release (*J*_{rel}^{NS}) generally decreases as *N* increases (Fig. 8*C*).

Figure 9 shows how steady-state probability density function, ρ_{ss}(*f*), and the release flux density, *v*_{rel}^{T}*f*(c̄_{jsr} − c̄_{ds})ρ_{ss}, depend on the number of RyRs per release site (*N*) when the total release rate *v*_{rel}^{T} is fixed (i.e., *MN* is a constant). For network SR [Ca^{2+}] of c_{nsr} = 950 μM (Fig. 9*A*), a larger number of channels per CaRU (*N*) decreases the “diffusion” term (channel gating fluctuations) in *Eq. 11* and both spark- and nonspark-mediated SR Ca^{2+} release. However, for a slightly larger value of c_{nsr} = 1,020 μM, larger *N* decreases nonspark-mediated release (*J*_{rel}^{NS}) while promoting robust sparks and increasing spark-mediated release *J*_{rel}^{S}. However, if the total release flux (*v*_{rel}^{T}) is proportional to *N* (as opposed to a constant), larger *N* results in higher release flux regardless of c_{nsr} because of high release flux rate (see appendix e).

Finally, Fig. 10 shows the steady-state spark *Score* for “clamped” c_{myo} and c_{nsr} and illustrates the interplay of bulk concentrations and Ca^{2+} sparks. For a given value of c_{myo}, the *Score* is a bell-shaped function of c_{nsr}, that is, there is a specific range of SR Ca^{2+} load that supports robust sparks. As observed in prior work (21), the range for robust sparks decreases as *N* is increased (Fig. 10, *B* and *C*). Most importantly, the solid black lines indicate the steady-state (unclamped) network SR [Ca^{2+}] (c_{nsr}) as a function of c_{myo} (cf. Fig. 4*B*). When c_{myo} is sufficiently elevated that further increase leads to decreased c_{nsr}, the SR Ca^{2+} load equilibrates to a value that maximizes the *Score*, that is, the steady-state c_{nsr} decreases (with increasing c_{myo}) just enough to maintain robust sparks. This intriguing and potentially significant result is also observed when the total release flux *v*_{rel}^{T} is proportional to *N* (not shown).

## DISCUSSION

#### Summary of main findings.

In this article, we present a novel local/global whole cell model of Ca^{2+} homeostasis based on a Langevin description of stochastic Ca^{2+} release that includes both spark-mediated and nonspark-mediated release dynamics. The Fokker-Planck equation associated with the Langevin formulation of stochastic Ca^{2+} release is coupled to balance equations for the bulk myoplasmic and network SR [Ca^{2+}]. With the use of this approximate representation of the collective dynamics of a large number of identical CaRUs, this whole cell modeling approach avoids Monte Carlo simulation of a large population of CaRUs and facilitates our study of Ca^{2+} homeostasis in permeabilized ventricular myocytes.

In permeabilized myocytes, the interplay between bulk myoplasmic [Ca^{2+}] (c_{myo}), and network SR [Ca^{2+}] (c_{nsr}) on SR Ca^{2+} release is complex, in spite of the fact that myoplasmic [Ca^{2+}] is under experimental control (i.e., c_{myo} is not a dynamic variable but a model parameter). Elevated c_{myo} promotes Ca^{2+} uptake into the network SR via the SERCA pump, and this may elevate c_{nsr}. On the other hand, high c_{myo} and high c_{nsr} both promote increased SR Ca^{2+} release and depletion of SR Ca^{2+}.

We use the Langevin and Fokker-Planck local/global whole cell model of a permeabilized ventricular myocyte to characterize the depletion of network SR [Ca^{2+}] (c_{nsr}) that occurs via both spark-mediated release and nonspark-mediated release, as well as dependency of SR Ca^{2+} load on myoplasmic [Ca^{2+}] (c_{myo}). In agreement with recent experimental work (3), we find that spark-mediated release increases exponentially as c_{myo} increases, while nonspark-mediated release increases linearly (Fig. 7).

The interplay among c_{myo}, c_{nsr}, and spark- and nonspark-mediated release in the local/global whole cell model generates several phenomena of Ca^{2+} homeostasis in permeabilized cells that are worth highlighting. For example, the model predicts the presence of two distinct stable steady states that lead to the same SR Ca^{2+} load, one with low myoplasmic [Ca^{2+}] and predominantly nonspark-mediated SR Ca^{2+} release and another with high myoplasmic [Ca^{2+}] and release that is primarily spark mediated (Fig. 5). Significantly, in our permeabilized ventricular myocyte model, for any clamped myoplasmic [Ca^{2+}] (c_{myo}) that is large enough to trigger spark-mediated release, the resulting spontaneous stochastic Ca^{2+} release tends to decrease the network SR Ca^{2+} load just enough to maintain robust Ca^{2+} sparks (Fig. 10). To our knowledge this potentially significant characteristic of Ca^{2+} homeostasis in permeabilized cells has not previously been identified.

#### Physiological significance.

Significant effort in recent years has been devoted to understanding the mechanisms influencing RyR regulation and SR Ca^{2+} release. Abnormal regulation of RyRs can lead to aberrant SR Ca^{2+} release that directly contributes to excitation-contraction coupling dysfunction (13, 14). Previous studies have shown RyR-mediated Ca^{2+} release was enhanced in myocytes from failing rabbit hearts (40), which increases the likelihood of Ca^{2+}-dependent arrhythmias (13). Recent experiments suggested that hidden RyR release contributes to the total release flux and influences Ca^{2+} homeostasis (3, 4, 40). In this article, we are particularly interested in how myoplasmic [Ca^{2+}] (c_{myo}) influences SR Ca^{2+} release via regulation of stochastic Ca^{2+} release mediated by CaRUs composed of clusters of RyRs. Our model shows that RyRs may produce both visible (spark-mediated) and invisible (nonspark-mediated) stochastic Ca^{2+} release. High c_{myo} increases both spark- and nonspark-mediated release by increasing the open probability of Ca^{2+}-activated RyRs. However, c_{myo} affects these pathways in two distinct and characteristic ways. Nonspark-mediated Ca^{2+} release increases linearly as a function of c_{myo}, while spark-mediated release increases exponentially with c_{myo}.

We investigated how the number of RyRs in each individual CaRU influences network SR Ca^{2+} depletion and stochastic Ca^{2+} release. When *v*_{rel}^{T} is fixed (single channel conductance inversely proportional to *N*), we found that a larger number of RyRs per CaRU results in a steeper release flux (primarily spark-mediated release) as a function of network SR [Ca^{2+}], when the SR is replete. However, when network SR [Ca^{2+}] is depleted, and the release flux is primarily nonspark mediated, increasing the number of RyRs per CaRU decreases the total release flux, due to reduced triggering of Ca^{2+} sparks (Figs. 8 and 9). When *v*_{rel}^{T} is proportional to *N* (fixed single channel conductance), SR Ca^{2+} decreases with increasing *N*, due to higher release rates (not shown).

Because recent studies have shown that the number of RyRs per CaRU is variable (1), we note that the local/global whole cell model presented here can be modified to account for CaRUs of different size by simultaneously solving multiple Fokker-Planck equations, each with a different value for *N*. Assuming *M* = ∑_{i}*M*_{i} CaRUs, with CaRUs of type *i* composed of *N*_{i} RyRs, the population densities ρ_{i} solve

where α_{i} = *v*_{i}^{+} − *v*^{−}, γ_{i} = (*v*_{i}^{+} + *v*_{i}^{−})/*N*, *v*_{i}^{+} = *k*^{+}(c̄_{ds}^{i})^{η}(1 − *f*) and *v*^{−} = *k*^{−}*f*. The stochastic Ca^{2+} release flux (*Eq. 27*) becomes

where ∫ρ_{i}d*f* = 1 and thus *M*^{−1}∑_{i}∫*M*_{i}ρ_{i}d*f* = 1. In these equations, c̄_{ds}^{i}(*f*) and c̄_{jsr}^{i}(*f*) are given by indexed versions of *Eqs. 32* and *33* where = *v*_{rel}^{0}*N*_{i}*f* and *v*_{rel}^{0} is analogous to the RyR unitary conductance. Writing *v*_{myo}^{i} and *v*_{nsr}^{i} as the domain time constants for a representative of the *i*th class of CaRU, χ_{myo}^{i} and χ_{nsr}^{i} are given by *Eqs. 34–37* upon replacement of *i* for T. The Fokker-Planck equations are coupled, because α_{i} is a function of c_{nsr} through c̄_{ds}^{i}, and dc_{nsr}/d*t* depends on the ρ_{i} through *v*_{nsr}^{T} (*Eq. 32*).

#### Comparison to other whole cell models.

A number of mathematical and computational whole cell models have been developed to understand Ca^{2+} homeostasis and the cardiac Ca^{2+} cycle. For example, computational models of excitation-contraction coupling in ventricular myocytes have been developed in which SR Ca^{2+} release depends directly on the average myoplasmic [Ca^{2+}] (25, 34). These “common pool” models (33) exhibit all-or-none triggered SR Ca^{2+} release, contrary to experiments showing that release is smoothly graded with changes in Ca^{2+} influx (5, 36). This discrepancy is a consequence of the “local control” mechanism of CICR. In ventricular myocytes, the cellular SR Ca^{2+} release flux is not a function of the spatially averaged intracellular [Ca^{2+}] but instead depends on thousands of different local Ca^{2+} concentrations fluctuating in response to stochastic openings and closings of RyRs located on the SR membrane. The picture is further complicated by dynamic changes in localized SR [Ca^{2+}] that are also spatially heterogeneous and thought to influence the gating of RyRs (31).

To overcome this problem, stochastic models that account for the heterogeneous dyadic subspace and junctional SR [Ca^{2+}] have been developed (19, 22, 39). Similar to the Langevin model that is the focus of this article, these local control models include a large number of CaRUs. In such models, RyR stochastic gating is typically described by a discrete-state Markov chain. This approach has recently been used to examine issues such as allosteric coupling between RyRs (39) and refractoriness of Ca^{2+} release after termination (28).

While Markov chain and Langevin models of CaRUs may lead to similar results (Fig. 1), the state space for Markov chain simulations is proportional to the number of CaRU states, a quantity that is exponential in the number of distinct RyR states. To see this, consider a CaRU composed *N* identical *K*-state channels (and thus *K*^{N} states). It is well-known that the number of distinguishable CaRU states is given by (*N* + *K* − 1)!/*N*!/(*K* − 1)! = [(*N* + *K* − 1)...(*N* + 1)]/(*K* − 1)!, a quantity that includes a term proportional to *N*^{K − 1} (the numerator has *K* − 1 terms) and is thus exponential in *K*. On the other hand, the run time for Langevin simulations is independent of the number of RyRs (*N* is a model parameter that scales the channel noise) and proportional to the number of RyR states *K* (the required number of SDEs). Similarly, the run time of the Langevin local/global model does not scale with *N*, and the model may be extended to include RyRs with more than two states (see below). Because the Langevin version of the local/global model that has been our focus involves only a single SDE (two-state RyR model), the probability density function for CaRU state is univariate. For this reason, the Fokker-Planck local/global whole cell model is extremely computationally efficient. Because a *K*-state RyR model leads to a Fokker-Planck equation with *K* − 1 independent variables (conservation of probability), the Langevin version of the local/global model is likely to be more straightforward than the Fokker-Planck version when *K* ≥ 3 (see *Eq. 41* below).

It is instructive to compare the local/global model presented here with our prior work. In Hartman et al. (22), we presented a similar minimal model of a permeabilized myocyte, in which bulk myoplasmic and network SR Ca^{2+} levels were coupled to a Markov chain CaRU model with *N* Ca^{2+}-activated RyRs per release site. The master equation in this case was a linear system of *N* + 1 ODEs. The Langevin and Fokker-Planck local/global models presented here are also distinct from prior work of Williams et al. (37, 38). In these studies, Ca^{2+} release dynamics were described by a set of coupled multivariate probability density functions (advection-reaction equations) for the dyadic subspace and junctional SR [Ca^{2+}], c_{ds} and c_{jsr}, conditioned on CaRU state. This population density method and the associated moment-based reductions (38) are limited by a state-space explosion that is exponential in *K*, while the computational efficiency of the Langevin local/global model is linear in *K*.

#### Limitations and extensions of the model.

In the Langevin model, we assume that the number of channels in each CaRUs is large enough that the fraction of RyRs in different states can be treated as a continuous variable. When the number of RyRs per CaRU is small, the error associated with the Langevin approximation to the Markov chain CaRU model may not be acceptable (15). In the local/global whole cell model presented here, the Langevin formulation was validated using a physiologically realistic numbers of RyRs per CaRU (see Figs. B1 and B2). The number of RyRs per CaRU required for the Langevin formulation to be highly accurate likely depends on the details of the RyR model used but is easily determined in any specific case.

In the derivation of the reduced local/global model, we assume that the dynamics of dyadic subspace [Ca^{2+}] and junction SR [Ca^{2+}] are fast compared with the gating of RyRs. However, slow translocation of junctional SR [Ca^{2+}] can be incorporated into the Langevin local/global whole cell model through the addition of an additional SDE (24). This extension might be important if the chosen RyR model includes luminal regulation, that is, transitions whose rate is a function of junctional SR [Ca^{2+}]. Accounting for slow junctional SR dynamics would increase the dimensionality of the probability density function (*Eq. 28*) used in the corresponding Fokker-Planck whole cell model.

Upgrading the Langevin formulation of the local/global whole cell model to accommodate more complex RyR models is straightforward. A *K*-state RyR model leads to a linear system of *K* SDEs,
(39)

where **f** = (*f*_{1}, *f*_{2}, . . ., *f*_{K}) and **ξ** = (ξ_{1}, ξ_{2}, . . ., ξ_{K}) are row vectors, *Q* = (*q*_{ij}) is the RyR model's transition matrix (the Markov chain's infinitesimal generator), the random term is mean zero ( = **0**) with two-time covariance matrix,
(40)

where Γ = (γ_{ij}), γ_{ij} = −(*q*_{ij}*f*_{i} + *q*_{ji}*f*_{j})/*N* for *i* ≠ *j* and γ_{ii} = −∑_{j≠i} γ_{ij} (the γ_{ii} are positive) (26). The corresponding Fokker-Planck equation for the *K*-state RyR is
(41)

where (**f**Q)_{i} is the *i*th element of the row vector **f**Q.

The Langevin and Fokker-Planck formulations of the local/global whole cell model presented here are not explicitly spatial. That is, a large population of CaRUs are assumed to influence one another indirectly via the spatially averaged bulk myoplasmic and network SR [Ca^{2+}] (the CaRUs are mean-field coupled). This form of the local/global model is not well-suited to investigate macrosparks and other explicitly spatial phenomena that might occur in permeabilized ventricular myocytes when myoplasmic [Ca^{2+}] is very high. By partitioning (discretizing) the bulk myoplasm and SR into regions that interact via buffered Ca^{2+} diffusion, the formalism would allow for propagation of intercellular Ca^{2+} waves, subcellular alternans, and so on. Such extensions of the Langevin local/global model approach would be straightforward and robust. Extending the Fokker-Planck local/global model in this way require a discretization sufficiently coarse that number of CaRUs per subcompartment remains large.

## GRANTS

The work was supported in part by National Science Foundation Grant DMS 1121606 (to G. D. Smith) and the Biomathematics Initiative at The College of William & Mary.

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the author(s).

## AUTHOR CONTRIBUTIONS

Author contributions: X.W., S.H.W., E.A.S., and G.D.S. conception and design of research; X.W. performed experiments; X.W. analyzed data; X.W., S.H.W., and G.D.S. interpreted results of experiments; X.W. prepared figures; X.W. drafted manuscript; X.W., S.H.W., Y.H., E.A.S., and G.D.S. edited and revised manuscript; X.W., S.H.W., Y.H., E.A.S., and G.D.S. approved final version of manuscript.

## Appendix A: DERIVATION OF THE LANGEVIN CaRU MODEL

Assuming a time interval Δ*t* is small enough so that at most one event occurs in the interval [*t*, *t* + Δ*t*], a CaRU composed of *N* two-state channels will undergo a *N*_{O} → *N*_{O} − 1, transition (Δ*N*_{O} = −1) with probability *p*^{−} = *k*^{−}*N*_{O}Δ*t*, and a *N*_{O} → *N*_{O} + 1 transition (Δ*N*_{O} = 1) with probability *p*^{+} = (*N* − *N*_{O})*k*^{+}c_{NO}^{η}Δ*t*, where c_{NO} is the local [Ca^{2+}] for *N*_{O} open channels (*Eq. 3*). By conservation, the probability that Δ*N*_{O} = 0 is *p*^{0} = 1 − *k*^{−}*N*_{O}Δ*t* − (*N* − *N*_{O})*k*^{+}c_{NO}^{η}Δ*t*. Conditioning on the current state, the expected infinitesimal increment of the number of open channels (Δ*N*_{O}) is thus

The deterministic part of the right hand side of *Eq. 5* is derived as the corresponding infinitesimal expected increment in *f*_{O} = *N*_{O}/*N*. Similarly, the infinitesimal variance of Δ*N*_{O} is

The function γ(*f*_{O}) that occurs in *Eq. 8* is derived from this quantity using E[Δ*f*_{O}^{2}|*f*_{O}(*t*)] = E[Δ*N*_{O}^{2}|*N*_{O}(*t*)]/*N*^{2}.

## Appendix B: COMPARISON OF MARKOV CHAIN AND LANGEVIN CaRU MODELS

Figure B1*A* compares the spark *Score* calculated via the Langevin (+ symbols) and the Markov chain (lines) description of a CaRU composed of two-state channels. The *Score* is a biphasic function of the coupling strength c_{*} (*Eq. 38*), with robust sparks occurring over a wider range of coupling strength when *N* = 20 vs. 60 (dashed and solid lines, respectively). The Langevin method agrees with the Markov chain result, but overestimates the *Score* slightly for *N* = 20 and small c_{*} (parameter regimes with few channel openings). Figure B1*B* shows that the *Score* calculated via the stationary distribution of the Markov chain and the Fokker-Planck equation are in agreement.

The Langevin method is also applicable to more complex single channel models. For example, consider a three-state RyR that is activated as well as inactivated by Ca^{2+},
(42)

where c is the local [Ca^{2+}], *k*_{a}^{+}c^{η}, *k*_{a}^{−}, *k*_{b}^{+}c^{η}, and *k*_{b}^{−} are transition rates with units of time^{−1}, *k*_{a}^{+} and *k*_{b}^{+} are association rate constants with units of conc^{−η}·time^{−1}, and the cooperativity of Ca^{2+} binding η is the same for activation and inactivation. The Langevin description of a CaRU composed of *N* three-state channels (*Eq. 42*) is given by *Eq. 39*, where the fraction of channels in each state, **f** = (*f*_{C}, *f*_{O}, *f*_{R}), is a row vector, *Q* is the transition rate matrix,

where the diagonal elements (◇) are such that each row sums to zero (∑_{j}*q*_{ij} = 0) and the local [Ca^{2+}] is c = c_{∞} + *f*_{O}c̄ where c̄ = *N*c_{*} (cf. *Eq. 3*). The rapidly varying forcing function in the Langevin equation (*Eq. 39*), **ξ**(*t*) = [ξ_{C}(*t*), ξ_{O}(*t*), ξ_{R}(*t*)] is mean zero ( = **0**) with two-time covariance = Γ(**f**)δ(*t* − *t*′) (*Eq. 40*). Here Γ = (γ_{ij}) is given by γ_{OC} = γ_{CO} = [*k*_{a}^{+}c^{η}*f*_{C} + *k*_{a}^{−}*f*_{O}]/*N*, γ_{OR} = γ_{RO} = [*k*_{b}^{+}c^{η}*f*_{O} + *k*_{b}^{−}*f*_{R}]/*N*, γ_{CR} = γ_{RC} and the diagonal entries are such that each row sums to zero.

Figure B2 plots *Score* vs. coupling strength (c_{*}) for this Langevin model of a CaRU composed of *N* three-state channels with Ca^{2+} inactivation. This may be compared with the result for a CaRU composed of *N* two-state channels with no inactivation (Fig. B1). Consistent with a previous computational study (21), Fig. B2 shows that Ca^{2}-dependent inactivation facilitates spark termination (i.e., CaRUs spark for a wider range of coupling strengths). Most importantly, the Langevin (+ symbols) and Markov chain (lines) simulations agree.

## Appendix C: LANGEVIN EQUATION BOUNDARY CONDITIONS

Because solutions of the Langevin CaRU model (*f*_{i}) represent the fraction of channels in state *i*, physical values are in the range 0 ≤ *f*_{i} ≤ 1 and, formally, the stochastic processes that solve the Langevin CaRU models (*Eqs. 5* and *39*) have this property. However, numerical integration via the Euler-Maruyama method (16) involves a finite time step; consequently, there is a small probability of crossing *f*_{i} = 0 or 1, thereby exiting the physical range.

In the context of stochastic ODE modeling of ion channel dynamics, several modifications of the Euler-Maruyama scheme are commonly used to address this numerical issue. These include rejection and projection methods as well as more sophisticated approaches such as equilibrium noise approximations (18) and reflected stochastic differential equations (reviewed in Ref. 7). Unfortunately, these methods yield solutions that may disagree with the corresponding Markov chains when *N* = 20–200 (11). In the context of Langevin CaRU models, a superior approach is to define auxiliary variables (observables) restricted to the physical range, i.e., = max[0, min(1, *f*_{i})], for evaluation of state-dependent rates, without projecting the stochastic trajectory *f*_{i} to the boundary. For example, the Euler-Maruyama scheme use to integrate *Eq. 5* is

where the Δ*B*^{m} are i.i.d. normal random variables with mean zero and variance 1/Δ*t*, . Because the deterministic flux is positive (α > 0) when *f* < 0 and negative (α < 0) when *f* > 1, no restriction is necessary for the factors 1 − *f* and *f* in ; in fact, we found that not doing so yields better agreement with the corresponding Markov chain simulations. Conversely, is used in the evaluation of the diffusive term to ensure . This method has similarities to the reflected stochastic differential equation technique discussed in Dangerfield et al. (2012).

## Appendix D: SPARK STATISTICS ANALYSIS VIA THE LANGEVIN DESCRIPTION OF THE LOCAL/GLOBAL WHOLE CELL MODEL

Figure D1 shows the mean steady-state spark amplitude (*A*), spark duration (*B*), and interevent intervals (*C*) as a function of c_{myo}, calculated via the Langevin version of the local/global whole cell model. The duration of the *i*th Ca^{2+} release event is the time elapsed between the first channel opening and last channel closing of each simulated spark, here defined as *f*_{O} crossing the threshold 1/*N* in the upward/downward direction. The amplitude of *i*th Ca^{2+} release event is the integrated area under *f*_{O}(*t*) during the event. The *i*th interevent interval is the length of time between the (*i* − 1)th and *i*th Ca^{2+} release events. Note that spark amplitude and spark duration are biphasic functions of c_{myo}, peaking at c_{myo} ≈ 0.25 μM, similar to the steady-state c_{nsr} and spark *Score* (Fig. 4, *B* and *C*).

## Appendix E: Ca^{2}^{+} RELEASE FLUX AND CaRU SIZE

Most of the parameter studies presented above assume that the total number of RyRs per cell is fixed. When the number of channels per CaRU (*N*) is varied, the number of CaRUs per cell (*M*) is changed so that *MN* is a constant (i.e., the total release flux rate *v*_{rel}^{T} is fixed). Alternatively, *M* may be fixed; in this case, *v*_{rel}^{T} is proportional to CaRU size (*N*). Figure E1 shows the total release flux (*J*_{rel}^{T}), spark-mediated release (*J*_{rel}^{S}), and nonspark-mediated release (*J*_{rel}^{NS}) when the number of channels per CaRU (*N*) are varied under this assumption (fixed single channel conductance). In this case, regardless of c_{nsr}, the total release flux and spark-mediated release are higher for larger *N*. Conversely, when *v*_{rel}^{T} is fixed (Fig. 8), the clamped network SR [Ca^{2+}] c_{nsr} determines whether CaRU size *N* increases or decreases the total release flux *J*_{rel}^{T}. Figure E2 shows release flux density increases with CaRU size when *v*_{rel}^{T} is proportional to *N* (cf. Fig. 9).

- Copyright © 2015 the American Physiological Society