This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.106.099861v1 92/7/2311    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Williams, G. S. B. Articles by Smith, G. D. Search for Related Content PubMed PubMed Citation Articles by Williams, G. S. B. Articles by Smith, G. D.

A Probability Density Approach to Modeling Local Control of Calcium-Induced Calcium Release in Cardiac Myocytes

George S. B. Williams ^*, 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; Medical Biotechnology Center and the Institute of Molecular Cardiology, University of Maryland Biotechnology Institute, Baltimore, Maryland; and Department of Pediatrics, New York University School of Medicine, New York, New York

Correspondence: Address reprint requests to Gregory D. Smith, Dept. of Applied Science, McGlothlin-Street Hall, Rm. 305, College of William and Mary, Williamsburg, VA 23187. E-mail: greg{at}as.wm.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS DISCUSSION APPENDIX A: DESCRIPTION OF... APPENDIX B: GENERALIZATION OF... APPENDIX C: DERIVATION OF... APPENDIX D: NUMERICAL SCHEME... ACKNOWLEDGEMENTS REFERENCES

 
We present a probability density approach to modeling localized Ca²⁺ influx via L-type Ca²⁺ channels and Ca²⁺-induced Ca²⁺ release mediated by clusters of ryanodine receptors during excitation-contraction coupling in cardiac myocytes. Coupled advection-reaction equations are derived relating the time-dependent probability density of subsarcolemmal subspace and junctional sarcoplasmic reticulum [Ca²⁺] conditioned on "Ca²⁺ release unit" state. When these equations are solved numerically using a high-resolution finite difference scheme and the resulting probability densities are coupled to ordinary differential equations for the bulk myoplasmic and sarcoplasmic reticulum [Ca²⁺], a realistic but minimal model of cardiac excitation-contraction coupling is produced. Modeling Ca²⁺ release unit activity using this probability density approach avoids the computationally demanding task of resolving spatial aspects of global Ca²⁺ signaling, while accurately representing heterogeneous local Ca²⁺ signals in a population of diadic subspaces and junctional sarcoplasmic reticulum depletion domains. The probability density approach is validated for a physiologically realistic number of Ca²⁺ release units and benchmarked for computational efficiency by comparison to traditional Monte Carlo simulations. In simulated voltage-clamp protocols, both the probability density and Monte Carlo approaches to modeling local control of excitation-contraction coupling produce high-gain Ca²⁺ release that is graded with changes in membrane potential, a phenomenon not exhibited by so-called "common pool" models. However, a probability density calculation can be significantly faster than the corresponding Monte Carlo simulation, especially when cellular parameters are such that diadic subspace [Ca²⁺] is in quasistatic equilibrium with junctional sarcoplasmic reticulum [Ca²⁺] and, consequently, univariate rather than multivariate probability densities may be employed.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS DISCUSSION APPENDIX A: DESCRIPTION OF... APPENDIX B: GENERALIZATION OF... APPENDIX C: DERIVATION OF... APPENDIX D: NUMERICAL SCHEME... ACKNOWLEDGEMENTS REFERENCES

 
The mechanical function of the heart depends on complex bidirectional interactions between electrical and calcium (Ca²⁺) signaling systems. Each time the heart beats, current flowing through the ion channels in the plasma membrane (sarcolemma) causes a characteristic change in membrane voltage known as an action potential (AP). Membrane depolarization during the AP causes L-type Ca²⁺ channels to open, and Ca²⁺ current through these channels causes the release of a larger amount of Ca²⁺ from the sarcoplasmic reticulum, a process known as Ca²⁺-induced Ca²⁺ release (CICR). This leads to a large, transient increase in [Ca²⁺] in each heart cell, and contraction occurs when these Ca²⁺ ions bind to myofilaments, a sequence of events known as excitation-contraction (EC) coupling. In addition, intracellular [Ca²⁺] feeds back upon and changes the cell's membrane potential through the Ca²⁺ dependence of several ion channels and membrane transporters.

Mathematical and computational modeling has proved to be an important tool for understanding cardiac electrophysiology and EC coupling. Computer simulations have been used to test hypotheses about heart cell function and predict underlying mechanisms (1–4). Most investigations have employed deterministic models that ignore molecular fluctuations and assume an isopotential cell, an approach that is valid for simulating current flowing through a large population of voltage-gated ion channels. Even though the individual channels open and close stochastically, each channel experiences the same voltage, so identical rate constants apply to each channel and only average behavior needs to be considered. However, this approach is not suitable for simulating CICR release during EC coupling because the overall release flux represents a collection of discrete events, known as Ca²⁺ sparks, evoked by local—rather than global—increases in Ca²⁺ concentration (5). That is, each spark reflects Ca²⁺ release from a cluster of Ca²⁺-regulated intracellular Ca²⁺ channels known as ryanodine receptors (RyRs) that is triggered by entry of Ca²⁺ through nearby L-type Ca²⁺ channels (6). Thus, different groups of RyRs experience different local Ca²⁺ concentrations and stochastically gate in a manner that depends on whether nearby sarcolemmal Ca²⁺ channels have recently been open or closed. One consequence of this "local control" (7) mechanism of cardiac CICR is that deterministic "common pool" models—whole cell models in which all RyR clusters in a myocyte experience the same [Ca²⁺]—fail to reproduce several important experimental observations. In particular, the high gain and positive feedback of common pool models ensures that Ca²⁺ is released in an all-or-none fashion (2,3,8–10) as opposed to being graded with the amount of Ca²⁺ influx, as observed in numerous experiments (6,11,12). Deterministic common pool models of cardiac CICR during EC coupling that have been able to reproduce graded release have done so in an ad hoc fashion (4,13–16).

Models of EC coupling are able to simulate graded Ca²⁺ release mechanistically by treating L-type Ca²⁺ channels and juxtaposed Ca²⁺ release sites as stochastic "Ca²⁺ release units" (CaRUs), each of which is associated with its own diadic subspace Ca²⁺ concentration. When activated spontaneously or through membrane depolarization these CaRUs may deplete Ca²⁺ stored in localized regions of junctional SR and, on a slower timescale, interact with one another via diffusion of Ca²⁺ within the network SR and bulk myoplasm. This approach, however, requires relatively large computational resources to perform Monte-Carlo simulations of stochastic Ca²⁺ release from a large population of CaRUs. Indeed, the number of simulated CaRUs is often reduced to unphysiological values in such models to obtain shorter run times (7,17–19).

Two recent deterministic models have used a minimal Ca²⁺ release unit formulation of interactions between L-type channels and RyR clusters to produce graded release (20,21). In these models ordinary differential equations for the fraction of Ca²⁺ release units in each of a small number of states are solved under the assumption that subspace [Ca²⁺] is an algebraic function of the bulk myoplasmic and network SR [Ca²⁺]. This function depends on Ca²⁺ release unit state and is determined by balancing the Ca²⁺ fluxes into and out of the diadic subspace. While the large number of Ca²⁺ release units in cardiac myocytes—estimated in the range of 10,000–20,000 via both structural (22) and functional (23) observations—does indeed suggest that it should be possible to produce deterministic local control models of EC coupling, the assumption that diadic subspace [Ca²⁺] is in quasistatic equilibrium with bulk myoplasmic and network SR Ca²⁺ may be overly restrictive. Indeed, this modeling approach is only valid when the dynamics of subspace [Ca²⁺] are very fast compared to stochastic Ca²⁺ release unit transition rates. Moreover, [Ca²⁺] in a particular subspace is likely to depend on the local "junctional" SR [Ca²⁺] rather than the bulk or network SR [Ca²⁺], especially if junctional SR depletion influences RyR gating, as suggested by both simulations (18) and recent experiments (24,25).

Here we present an alternative deterministic formalism for modeling local control of CICR during cardiac EC coupling that captures the collective behavior of a large population of Ca²⁺ release units without this restrictive assumption. We utilize the fact that the number of Ca²⁺ release units is large (similar to Hinch (20) and Greenstein et al. (21)), but we do not assume a simple algebraic relationship between the local diadic subspace [Ca²⁺] associated with each Ca²⁺ release unit and the bulk Ca²⁺ concentrations. Instead, we define a set of multivariate continuous probability density functions for the diadic subspace and junctional SR [Ca²⁺] conditioned on CaRU state (26–28). As described below, these probability density functions solve a system of advection-reaction equations that are derived from the stochastic ordinary differential equations used in Monte Carlo simulations of local control. These equations are solved numerically using a high-resolution finite difference scheme while coupled to ordinary differential equations for the bulk myoplasmic and network SR [Ca²⁺]. This produces a minimal model of cardiac EC coupling that avoids computationally demanding Monte Carlo simulation while accurately representing heterogeneous local Ca²⁺ signals; in particular, the statistical recruitment of CaRUs and the dynamics of junctional SR depletion, spark termination, and junctional SR refilling.

Some of these results have previously appeared in abstract form (29).

	MODEL FORMULATION

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS DISCUSSION APPENDIX A: DESCRIPTION OF... APPENDIX B: GENERALIZATION OF... APPENDIX C: DERIVATION OF... APPENDIX D: NUMERICAL SCHEME... ACKNOWLEDGEMENTS REFERENCES

 
The minimal whole cell model of cardiac EC coupling that is the focus of this article can be formulated as a traditional Monte Carlo calculation in which heterogeneous local Ca²⁺ signals associated with a large number of CaRUs are simulated. In this Monte Carlo formulation, a diadic subspace and junctional SR compartment is associated with each CaRU and the [Ca²⁺] in these compartments is found by solving a large number of ordinary differential equations. Alternatively, these heterogeneous local Ca²⁺ signals can be simulated using a novel probability density approach that represents the distribution of diadic subspace and junctional SR Ca²⁺ concentrations with a system of partial differential equations (see below). Because many of the equations and parameters of the whole cell model of EC coupling are identical in the two formulations, we begin by presenting the Monte Carlo formulation.

Whole cell model of EC coupling: Monte Carlo formulation
Fig. 1 shows a diagram of the components and fluxes of the model of local Ca²⁺ signaling and CaRU activity during cardiac EC coupling that is the focus of this article. As illustrated in Fig. 1 A, each Ca²⁺ release unit includes two restricted compartments (the diadic subspace and junctional SR) with [Ca²⁺] denoted by , respectively, where the superscripted n is an index over a total number of Ca²⁺ release units (denoted by N). Each Ca²⁺ release unit includes an L-type Ca²⁺ channel dihydropiridine receptor (DHPR) and a minimal representation of a cluster of RyRs that is either fully closed or fully open. The fluxes indicate Ca²⁺ entry into a subspace via the DHPR or RyR cluster, respectively. Diffusion of Ca²⁺ between the nth diadic subspace and bulk myoplasm (c_myo) is indicated by . Similarly, indicates diffusion between the network SR (c_nsr) and junctional SR compartment associated with the nth Ca²⁺ release unit.

View larger version (24K): [in this window] [in a new window]   FIGURE 1  Diagrams of model components and fluxes. (A) Each Ca2+ release unit consists of two restricted compartments (the diadic subspace and junctional SR with [Ca2+] denoted by cds and cjsr, respectively), a two-state L-type Ca2+ channel (DHPR), and a two-state Ca2+ release site (a RyR "megachannel" (18)). The t-tubular [Ca2+] is denoted by cext and the fluxes , , , , Jin, Jncx, Jserca, and Jleak are described in the text. (B) The bulk myoplasm (cmyo) and network SR (cnsr) Ca2+ concentrations in the model are coupled via to a large number of Ca2+ release units (for clarity only four are shown), each with different diadic subspace () and junctional SR () Ca2+ concentration.

A complete description of CICR would include stochastic gating of roughly N = 20,000 CaRUs, each of which would contain multiple L-type Ca²⁺ channels (1–10) (30) and RyRs (30–300) (31), with each individual channel described by a Markov chain that consists of two to several tens of states. However, previous Monte Carlo simulations of EC coupling focusing on local control have often used Markov models of reduced complexity (7,18,20). Because such minimal models capture the essential characteristics of EC coupling gain and gradedness in simulated whole cell voltage clamp protocols, this level of resolution will suffice for our main purpose, which is to introduce the probability density approach as an alternative to Monte Carlo simulation.

A minimal four-state Ca²⁺ release unit model
Previous modeling studies indicate that the gating of the cluster of RyRs associated with each CaRU is all-or-none (7,17,18) and this suggests the following minimal two-state model of an RyR "megachannel",

,

(1)

where the Ca²⁺ activation of the cluster of RyRs is a sigmoidal function of the diadic subspace [Ca²⁺] (18),

and the influence of junctional SR [Ca²⁺] on RyR gating is included by making the half-maximal activation of the RyR megachannel (K_ryr) a decreasing function of ,

so that depletion of the junctional SR will render CaRUs refractory to activation after release terminates (18).

Similarly, to illustrate and validate the probability density approach it is sufficient to consider a two-state model of the L-type Ca²⁺ channel (DHPR),

,

(2)

with a voltage-dependent activation rate given by (4)

and constant deactivation rate that sets the mean open time (0.2 ms) and maximum open probability (0.1) of the channel. Although this two-state DHPR model ignores voltage- and Ca²⁺-dependent inactivation of L-type Ca²⁺ channels, these processes do not significantly influence the triggering of CICR during the whole-cell voltage clamp protocols that are the focus of this article (cf. Hinch (20)).

When the kinetic schemes of the RyR megachannel and DHPR (Eqs. 1 and 2) are combined we obtain the following minimal four-state model of a Ca²⁺ release unit,

(3)

where the horizontal transitions represent RyR opening and closing whereas vertical transitions represent DHPR gating.

Concentration balance equations
In the Monte Carlo formulation of the minimal whole cell model of EC coupling there are 2 + 2N ordinary differential equations representing Ca²⁺ concentration balance for the bulk myoplasm, network SR, N diadic subspaces, and N junctional SRs. Consistent with Fig. 1 these equations are

(4)

(5)

(6)

(7)

where 1 n N in Eqs. 5 and 6 and the total efflux and refill fluxes occurring in Eqs. 4 and 7 include a contribution from each CaRU and thus are given by and . Similarly, the total (trigger) flux via DHPR channels and the total release flux via RyR megachannels throughout the whole cell model are given by

(8)

The effective volume ratios _nsr, _ds, and _jsr in Eqs. 5–7 are defined with respect to the physical volume (V_myo) and include a constant-fraction Ca²⁺ buffer capacity for the myoplasm (β_myo). For example, the effective volume ratio associated with the network SR is

with effective volumes defined by and . Because each individual diadic subspace is assumed to have the same physical volume (V_ds) and buffering capacity (β_ds), the effective volume ratio that occurs in Eq. 5 is

(9)

where the second expression defines _ds in terms of the total physical volume of all the diadic subspaces in aggregate (). Similar assumptions and equations apply for the junctional SR so that the definition of _jsr follows Eq. 9.

We also define an overall myoplasmic [Ca²⁺] that includes contributions from the bulk myoplasm and each of the N diadic spaces (scaled by their effective volumes),

(10)

where the second equality uses natural definitions for the total effective diadic subspace volume, , and the average diadic subspace [Ca²⁺],

(11)

Similarly, the overall SR [Ca²⁺] involves both the junctional and network SR,

(12)

where , , and the average junctional SR [Ca²⁺] is defined as .

Description of fluxes
The trigger Ca²⁺ flux into each of the N diadic spaces through DHPR channels ( in Eq. 5) is given by

(13)

where A_m = C_mβ_myo/V_myo. The inward Ca²⁺ current () is given by

(14)

where V = RT/zF, is the total (whole cell) permeability of the L-type Ca²⁺ channels, and is a random variable that is 0 when the L-type Ca²⁺ channel associated with the nth CaRU is closed and 1 when this channel is open (Eqs. 2 and 3).

Similarly, the flux through the RyR megachannel associated with the nth CaRU () is given by

(15)

where = 0 or 1 when the release site is closed or open, respectively (Eqs. 1 and 3). Diffusion from each subspace into the bulk myoplasm is given by

(16)

and, similarly, diffusion from the network SR to each junctional SR compartment is given by

(17)

The remaining four fluxes that appear in Eqs. 4–6 include J_in (background Ca²⁺ influx), J_ncx (Na⁺-Ca²⁺ exchange), J_serca (SR Ca²⁺-ATPases), and J_leak (the network SR leak). The functional form of these four fluxes that directly influence the bulk myoplasmic [Ca²⁺] follows previous work (3,32,33) (see Appendix A).

Whole cell model of EC coupling: probability density formulation
The probability density approach to modeling local Ca²⁺ signaling and CaRU activity during cardiac EC coupling is an alternative to Monte Carlo simulation that is valid when the number of Ca²⁺ release units is large. We begin by defining continuous multivariate probability density functions for the diadic subspace () and junctional SR () Ca²⁺ concentrations jointly distributed with the state of the Ca²⁺ release unit () (34,35,26), that is,

(18)

where the index runs over the four Ca²⁺ release unit states (see Eq. 3) and the tildes on , , and indicate random quantities. If the meaning of Eq. 18 is not obvious, it may be helpful to imagine performing a Monte Carlo simulation as described in the previous section with a very large number of CaRUs. At any time t one could randomly sample one CaRU from this population to produce an instance of the random variables , , and , corresponding to the current state of the sampled L-type channel and RyR cluster and the diadic subspace and junctional SR [Ca²⁺] associated with this CaRU. The quantity ⁱ(c_ds, c_jsr, t) defined in Eq. 18 simply indicates the probability with which you would find this sampled CaRU in state i with diadic subspace [Ca²⁺] in the range [c_ds, c_ds + dc_ds] and junctional SR [Ca²⁺] in the range [c_jsr, c_jsr + dc_jsr] provided the total number of CaRUs is very large.

For the multivariate probability densities defined by Eq. 18 to be consistent with the dynamics of the Monte Carlo model of cardiac EC coupling described in the previous section, they must satisfy the following system of advection-reaction equations (26–28),

(19)

(20)

(21)

(22)

where the advection rates are functions of c_ds and c_jsr that can be read off the ordinary differential equations for the evolution the diadic subspace and junctional SR [Ca²⁺]. Consistent with Eqs. 5 and 6 we have

(23)

(24)

where indicates whether or not the L-type Ca²⁺ channel is open (, ) and, similarly, indicates whether or not the RyR channel cluster is open (, ). Eqs. 23 and 24 include four fluxes that may influence the diadic subspace and junctional SR [Ca²⁺] and consistent with Eqs. 13–17 these are given by

(25)

(26)

(27)

(28)

The advection terms in Eqs. 19–22 involving partial derivatives with respect to c_ds and c_jsr correspond to the deterministic dynamics of diadic subspace and junctional SR Ca²⁺ that depend on Ca²⁺ release unit state via and (Eqs. 5 and 6). Conversely, the reaction terms in Eqs. 19 and 22 correspond to the stochastic gating of the four-state Ca²⁺ release unit model whose transition rates are presented above (Eqs. 1–3). That is, Ca²⁺ release unit state changes move probability from one joint probability density to another in a manner that may [] or may not [, , and ] depend on the diadic subspace and junctional SR [Ca²⁺].

It is important to note that the functional form of the fluxes and occurring in Eqs. 23 and 24 involve the bulk myoplasmic and network SR Ca²⁺ concentrations (c_myo(t) and c_nsr(t) in Eqs. 26 and 27). These bulk Ca²⁺ concentrations satisfy ordinary differential equations (ODEs) that are similar in form to the concentration balance equations used in the Monte Carlo approach (Eqs. 4 and 7),

(29)

(30)

where J_leak, J_ncx, J_serca, and J_in are defined as in the Monte Carlo approach (see Appendix A), but J_efflux* and J_refill* are functionals of the probability densities [ⁱ(c_ds, c_jsr, t)] governed by Eqs. 19–22, that is,

(31)

(32)

where is the total probability distribution of the diadic subspace and junctional SR [Ca²⁺] irrespective of the state of a randomly sampled CaRU, and the double integrals account for all possible values of diadic and junctional SR [Ca²⁺].

Summary of model formulation
The probability density and Monte Carlo formulations of the minimal model of EC coupling presented above have much in common. For example, the dynamics of the bulk myoplasmic and network SR [Ca²⁺] take similar forms (compare Eqs. 29 and 30 to Eqs. 4 and 7). However, the two approaches differ fundamentally in how the heterogeneous localized Ca²⁺ concentrations associated with a large number of Ca²⁺ release units are represented. In the traditional Monte Carlo simulation, 2N ordinary differential equations are solved to determine the dynamics of [Ca²⁺] in the diadic subspace and junctional SR compartments associated with N Ca²⁺ release units (Eqs. 5 and 6). In the probability density formulation, time-dependent multivariate probability densities for the diadic subspace and junctional SR [Ca²⁺] jointly distributed with CaRU state are updated by solving four coupled advection-reaction equations (Eqs. 19–22), one for each state of the chosen CaRU model (Eq. 3). Further details of the probability density approach are presented in Appendices B–D.

	RESULTS

TOP ABSTRACT INTRODUCTION MODEL FORMULATION RESULTS DISCUSSION APPENDIX A: DESCRIPTION OF... APPENDIX B: GENERALIZATION OF... APPENDIX C: DERIVATION OF... APPENDIX D: NUMERICAL SCHEME... ACKNOWLEDGEMENTS REFERENCES

 
In the following sections, traditional Monte Carlo simulations of voltage-clamp protocols using the minimal whole cell model of EC coupling presented above are shown to produce high-gain Ca²⁺ release that is graded with changes in membrane potential, a phenomenon not exhibited by so-called "common pool" models of excitation-contraction coupling. Analysis of these Monte Carlo results suggests a simplification of the advection-reaction equations that form the basis of the probability density approach. This reduced probability density formulation is subsequently validated against, and benchmarked for computational efficiency by comparison to, traditional Monte Carlo simulations.

Representative Monte Carlo simulations
Fig. 2 A shows representative Monte Carlo simulations of the minimal whole cell model of EC coupling presented above (Eqs. 1–17 and Appendix A). In this simulated voltage-clamp protocol, the holding potential of –80 mV is followed by a 20-ms duration test potential to –30, –20, and –10 mV (dotted, dot-dashed, and solid lines, respectively). Because these simulations involve a large but finite number of Ca²⁺ release units (N = 5000), the resulting Ca²⁺ influx through L-type Ca²⁺ channels (), elevation in the average diadic subspace concentration (), and the induced Ca²⁺ release flux () are erratic functions of time. As expected, the test potential of –10 mV leads to greater Ca²⁺ influx, higher diadic subspace [Ca²⁺], and more Ca²⁺ release than the test potentials of –30 and –20 mV. When the test potential is –10 mV a 30x "gain" is observed, here defined as the ratio where the overbar indicates an average over the duration of the pulse. Importantly, Ca²⁺ release exhibited by this Monte Carlo model is graded with changes in membrane potential (compare traces) and depolarization duration (not shown), phenomena that are not exhibited by common pool models of excitation-contraction coupling.

View larger version (12K): [in this window] [in a new window]   FIGURE 2  (A) Monte Carlo simulation of the whole cell model exhibits graded release during step depolarization from a holding potential of –80 mV to –30, –20, and –10 mV (dotted, dot-dashed, and solid lines, respectively). From top to bottom: command voltage, average diadic subspace [Ca2+] (, Eq. 11), total Ca2+ flux via L-type PM Ca2+ channels (, Eqs. 8, 13, and 14), and total Ca2+-induced Ca2+ release flux (, Eqs. 8 and 15). The simulation used N = 5000 Ca2+ release units. (B) Monte Carlo simulations similar to panel A except that the step potential is –10 (solid lines) and +10 mV (dotted lines), respectively. Here and below parameters are as in Tables 1–3.

The solid lines of Fig. 3 show [Ca²⁺] in the bulk myoplasm (c_myo) and network SR (c_nsr) during and after the –10 mV voltage pulse (note change in timescale). Approximately 400 ms is required for the bulk myoplasm and network SR concentrations to return to resting levels. Note that although the voltage pulse ends at t = 30 ms, the bulk myoplasmic [Ca²⁺] continues to increase for 20 ms. Similarly, the network SR [Ca²⁺] concentration continues to decrease until t = 80 ms.

View larger version (7K): [in this window] [in a new window]   FIGURE 3  Solid lines show the dynamics of bulk myoplasmic (cmyo) and network SR (cnsr) [Ca2+] in the whole cell voltage clamp protocol of Fig. 2 with step potential of –10 mV (note longer timescale). Dashed lines show the overall myoplasmic (, Eq. 10) and network SR (, Eq. 12) [Ca2+] that include contributions from diadic subspaces and junctional SR, respectively. Note that is only slightly greater than cmyo and the two traces are not distinguishable.

Dynamics of a representative Ca²⁺ release unit
Fig. 4 shows the dynamics of an individual Ca²⁺ release unit from the Monte Carlo simulations above (test potential of –10 mV, solid line of Fig. 2). Fig. 4 A shows the state of this representative Ca²⁺ release unit and the associated diadic subspace and junctional SR Ca²⁺ concentrations. When the DHPR initially opens (transition from state in Eq. 3) an influx of trigger Ca²⁺ leads to 7 µM increase in diadic subspace [Ca²⁺] and causes the RyR cluster to open ( transition). The resulting Ca²⁺-induced Ca²⁺ release quickly drives the diadic subspace [Ca²⁺] to 150 µM but over the next 10 ms the resulting decrease in junctional SR [Ca²⁺] leads to decreasing diadic subspace [Ca²⁺]. Note that junctional SR depletion is nearly complete in Fig. 4 before the transition that ends Ca²⁺ release; however, this example is not representative in this regard because most sparks terminate via stochastic attrition whereas depletion is only partial. Superimposed on the gradual decrease in diadic subspace [Ca²⁺] are square pulses of increased [Ca²⁺] (±7 µM) due to the stochastic openings of the L-type Ca²⁺ channel associated with this CaRU ( transitions).

View larger version (12K): [in this window] [in a new window]   FIGURE 4  (A) Dynamics of the diadic subspace () and junctional SR () Ca2+ concentrations associated with a single Ca2+ release unit during the voltage clamp protocol of Figs. 2 and 3. (B) The dynamics of these local Ca2+ concentrations in the (cds,cjsr)-plane. Trajectory color indicates CaRU state: both the L-type channel and the RyR cluster closed (, black); L-type channel open and RyR cluster closed (, green); L-type channel closed and RyR cluster open (, blue); both the L-type channel and the RyR cluster open (, red). Colored dashed lines correspond to estimates of diadic subspace [Ca2+] given by Eq. 33.

(33)

where and depend on Ca²⁺ release unit state and and are functions of plasma membrane voltage defined by with as in Eq. 28.

Fig. 4 B replots the dynamics of the diadic subspace and junctional SR [Ca²⁺] shown in Fig. 4 A in the (c_ds, c_jsr)-plane. The black arrows indicate the direction of the trajectories and color of the solid lines indicates CaRU state (, black; , green; , red; , blue). The diagonal trajectory is one consequence of diadic subspace [Ca²⁺] being "slaved" to junctional SR [Ca²⁺] as the junctional SR depletes. The four colored dotted lines correspond to the four functional relationships between and given by Eq. 33 (one for each CaRU state). The dynamics of diadic subspace [Ca²⁺] (solid lines) are well approximated by these dotted lines (save for short time intervals immediately following CaRU state transitions), demonstrating the validity of the quasistatic approximation leading to Eq. 33.

Dynamics of the population of Ca²⁺ release units
Fig. 4 shows the dynamics of the diadic subspace and junctional SR [Ca²⁺] associated with a single Ca²⁺ release unit during a voltage clamp step (Figs. 2 and 3). Conversely, Fig. 5 presents the state of each of the 5000 CaRUs at a particular moment in time (t = 30 ms, halfway through the test potential of –10 mV). To interpret this figure, it is important to understand that the four central panels of Fig. 5 correspond to the four CaRU states and are arranged in a manner corresponding to the transition state diagram of Eq. 3. At this moment during the simulation, 5% of the Ca²⁺ release units have open L-type channels () whereas 30% have an open RyR cluster (). Note that for each of the four subpopulations of CaRUs there is a linear relationship between c_ds and c_jsr, that is, the open circles tend to be arranged in lines, the position of which depends on CaRU state (and the slope of which depends on whether or not the RyR cluster is open). Thus, Fig. 5 demonstrates that across the entire population of Ca²⁺ release units, the observed diadic subspace [Ca²⁺] is well approximated by the quasistatic approximation given by Eq. 33.

View larger version (14K): [in this window] [in a new window]   FIGURE 5  The open circles are a snapshot at t = 30 ms of the diadic subspace () and junctional SR () Ca2+ concentrations in the Monte Carlo simulation of Fig. 2. Each of the four central panels corresponds to a particular Ca2+ release unit state and size of each subpopulation at this moment is indicated by through . The horizontally (vertically) oriented histograms give the marginal distribution of diadic subspace (junctional SR) [Ca2+] conditioned on CaRU state. Histograms are scaled for clarity and in some cases also truncated (asterisks).

A univariate probability density formulation for junctional SR [Ca²⁺]
It is important to note that the Monte Carlo simulations presented in Fig. 5 are only a snapshot of the population of 5000 Ca²⁺ release units. As the simulation progresses, imagine the open circles moving around in these four (c_ds, c_jsr)-planes consistent with Eqs. 5 and 6 with occasional jumps from one plane to another when a CaRU changes state. These four planes are analogous to the four time-dependent joint probability densities that form the basis of the probability density approach presented above (Eq. 18).

The observation that the diadic subspace [Ca²⁺] is well approximated by Eq. 33 across the entire population of Ca²⁺ release units (Fig. 5) suggests that the multivariate joint probability density functions defined in Eq. 18 will be well approximated by

(34)

where is a function of CaRU state and the junctional SR, bulk myoplasmic, and network SR [Ca²⁺] analogous to Eq. 33,

(35)

where , , , and are as defined in the previous section. The univariate probability density (c_jsr, t) that appears in Eq. 34 is the marginal density of the junctional SR [Ca²⁺] jointly distributed with CaRU state defined by

(36)

That is, when the observed form of the joint multivariate probability densities (Eq. 34) is integrated with respect to diadic subspace [Ca²⁺] we obtain

(37)

where the last equality uses the unit mass of the function, .

As shown in Appendix C, the observed form of the multivariate probability densities (Eq. 34) and the definition of the marginal density (first equality in Eq. 37) can be used to reduce Eqs. 19–22 into a univariate version of the probability density formulation that focuses on the dynamics of the marginal densities for the junctional SR [Ca²⁺] jointly distributed with CaRU state [(c_jsr, t)]. The resulting advection-reaction equations are (26–28),

(38)

(39)

(40)

(41)

where the advection rates , , , and are given by Eq. 24 with the substitution of for c_ds, that is,

(42)

(43)

where is the function of c_myo(t), c_jsr, and CaRU state (i) given by Eq. 35.

In this univariate probability density formulation, the bulk myoplasmic and network SR [Ca²⁺] are still given by Eqs. 29 and 30, but J_efflux* and J_refill* are now functionals of the joint marginal probability densities [(c_jsr, t)],

(44)

(45)

Comparison of probability density and Monte Carlo results
The four histograms presented in Fig. 6, A–D, show the marginal distributions of junctional SR [Ca²⁺] observed in Fig. 5 on identical scales. When presented in this fashion it becomes apparent that at t = 30 ms only a small fraction (5%) of the Ca²⁺ release units have open L-type Ca²⁺ channels (states and ), while 30% contain open RyR clusters (). Note that in Fig. 6 A the histogram bin representing Ca²⁺ release units with closed L-type Ca²⁺ channel, closed RyR cluster, and replete junctional SR is truncated; in fact, 80% of CaRUs in state have . With this understanding, a comparison of Fig. 6, A and B, shows that CaRUs with open RyR clusters are more likely to be depleted than CaRUs with closed RyR clusters, but CaRUs with closed RyR clusters are not necessarily replete, because recovery of junctional SR [Ca²⁺] is not complete until 400 ms after RyR closure (cf. Fig. 3).

View larger version (11K): [in this window] [in a new window]   FIGURE 6  Histograms of the junctional SR Ca2+ concentrations () at t = 30 ms in the Monte Carlo simulation of Figs. 2–5 jointly distributed with CaRU state. These histograms are plotted on the same scale, but one is truncated for clarity. For comparison, the solid lines show the four joint probability densities , , , and for junctional SR [Ca2+] (Eq. 34) calculated via numerical solution of Eqs. 29, 30, and 38–45. The probability density calculation of the fraction of subunits in each of the four states is denoted by i (Eq. 46).

(46)

in the probability density calculation is consistent with the Monte Carlo simulation Fig. 5, for example, in Fig. 6 A and this corresponds to in Fig. 5 A.

While Fig. 6 shows the four marginal probability densities [(c_jsr, t)] for the junctional SR [Ca²⁺] jointly distributed with CaRU state at a particular moment in time, Fig. 7 A shows the total probability density

(47)

evolving over time. Initially the mass of this probability density is concentrated at c_jsr 1000 µM (a in Fig. 7). During the 20-ms voltage pulse, a significant fraction of the probability density (65%) moves to junctional SR Ca²⁺ concentrations that are more than half-depleted (b), while 35% remains above 500 µM. Interestingly, the probability density remains bimodal for 200 ms after the voltage pulse ends (c and d). During this time, the probability mass that corresponds to depleted junctional SR (c) gradually moves to higher values of c_jsr as these junctional SR compartments are refilled via Ca²⁺ transport from the network SR. At the same time, the probability mass that corresponds to replete junctional SR compartments (d) follows the network SR [Ca²⁺] that decreases from t = 30–100 ms and increases again when t > 100 ms (recall the solid line in Fig. 3). Perhaps most importantly, Fig. 7 shows that the shape and temporal evolution of the distributions that form the basis of the probability density approach can be quite complicated.

View larger version (18K): [in this window] [in a new window]   FIGURE 7  Waterfall plot (A) and snapshots (B) of the time evolution of the total probability density for the junctional SR [Ca2+] (T(cjsr, t) given by Eq. 47) calculated via numerical solution of Eqs. 29, 30, and 38–45. The solid black lines show the 20-ms voltage step to –10 mV. See text for description of a–d.

Fig. 8 clarifies this point by showing how the total release flux (, open squares) observed in Monte Carlo simulation converges to the probability density result (solid lines) as the number of Ca²⁺ release units is increased from N = 50–20,000. Each panel shows a representative Monte Carlo simulation with voltage step to –10 mV (solid gray line) as well as the mean and standard deviation of 10 trials (open squares and error bars). As expected, the fluctuations in the total release flux decrease in magnitude as the number of CaRUs used in the Monte Carlo calculation increases. Similarly, Fig. 9 shows histograms of the junctional SR [Ca²⁺] (irrespective of CaRU state) at t = 30 ms in Monte Carlo simulations performed with a greater or lesser number of CaRUs. Notice that the probability density function ^T(c_jsr, t) (Eq. 47) accurately represents the distribution of junctional SR [Ca²⁺] so long as the number of CaRUs is 5000 or more. Indeed, in both Figs. 8 and 9 the Monte Carlo simulations are converging to the probability density result well before the Monte Carlo calculations include a physiological number of Ca²⁺ release units (N = 20,000). This indicates that the probability density approach to modeling local Ca²⁺ signaling and Ca²⁺ release unit activity in cardiac myocytes is a viable alternative to Monte Carlo simulation.

View larger version (17K): [in this window] [in a new window]   FIGURE 8  Total Ca2+ release flux () in Monte Carlo simulations utilizing increasing numbers of Ca2+ release units (N = 50, 500, 5000, and 20,000, respectively). Each panel shows a representative Monte Carlo simulation (solid gray line) and the mean and standard deviation of 10 trials (open squares and error bars). The solid lines show the corresponding probability density result (same in each panel).

View larger version (16K): [in this window] [in a new window]   FIGURE 9  Histograms of junctional SR [Ca2+] () at t = 30 ms in the Monte Carlo simulations similar to Fig. 5 but with increasing numbers of Ca2+ release units (N = 50, 500, 5000, and 20,000, respectively) One bin representing 57% probability of a replete junctional SR is truncated for clarity (asterisk). The solid lines show the probability density calculation of T(cjsr, t) (Eq. 47), the distribution of the total probability density for the junctional SR [Ca2+] (same in each panel).