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

Article Outline

Abstract
Introduction
Model formulation
Results
Discussion
Acknowledgments
Appendix A
Appendix B
Appendix C
Appendix D
References
Publication Information

Article Information

PubMed

Moment Closure for Local Control Models of Calcium-Induced C...

Termination of Cardiac Ca Sparks: An Investigative Mathemati...

Modeling Gain and Gradedness of Ca Release in the Functional...

…more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 92, Issue 7, 2311-2328, 1 April 2007

doi:10.1529/biophysj.106.099861

Biophysical Theory and Modeling

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
^§ Department of Pediatrics, New York University School of Medicine, New York, New York

Address reprint requests to Gregory D. Smith, Dept. of Applied Science, McGlothlin-Street Hall, Rm. 305, College of William and Mary, Williamsburg, VA 23187.

Abstract

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

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,2,3,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,9,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,14,15,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,18,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,27,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

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 Figure 1A, 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.

Display large version of this figure

Figure 1

Diagrams of model components and fluxes. (A) Each Ca²⁺ release unit consists of two restricted compartments (the diadic subspace and junctional SR with [Ca²⁺] denoted by c_ds and c_jsr, respectively), a two-state L-type Ca²⁺ channel (DHPR), and a two-state Ca²⁺ release site (a RyR “megachannel” 18). The t-tubular [Ca²⁺] is denoted by c_ext and the fluxes

, J_in, J_ncx, J_serca, and J_leak are described in the text. (B) The bulk myoplasm (c_myo) and network SR (c_nsr) Ca²⁺ concentrations in the model are coupled via

to a large number of Ca²⁺ release units (for clarity only four are shown), each with different diadic subspace (

) and junctional SR (

) Ca²⁺ concentration.

Figure 1B illustrates how the bulk myoplasm and network SR Ca²⁺ concentrations in the model are coupled via the diffusion fluxes () to a large number of Ca²⁺ release units (for clarity only four are shown). Importantly, each of the N Ca²⁺ release units may have a different diadic subspace () and junctional SR () Ca²⁺ concentration. Four additional fluxes directly influence the bulk myoplasm: a background Ca²⁺ influx denoted by J_in, extrusion of Ca²⁺ via the Na⁺-Ca²⁺ exchanger (J_ncx), SR Ca²⁺-ATPase (SERCA) pumps (J_serca) that resequester Ca²⁺ into the network SR, and a passive leak out of the network SR to the bulk myoplasm (J_leak).

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.2ms) 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)) 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 the total efflux and refill fluxes occurring in Eqs. (4) 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) 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)).

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

(15)

where

or 1 when the release site is closed or open, respectively (Eqs. (1)). 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) 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,27,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) 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) include four fluxes that may influence the diadic subspace and junctional SR [Ca²⁺] and consistent with Eqs. (13) these are given by

(25)

(26)

(27)

(28)

The advection terms in Eqs. (19) 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)). Conversely, the reaction terms in Eqs. (19) correspond to the stochastic gating of the four-state Ca²⁺ release unit model whose transition rates are presented above (Eqs. (1)). 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) involve the bulk myoplasmic and network SR Ca²⁺ concentrations (c_myo(t) and c_nsr(t) in Eqs. (26)). 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)),

(29)

(30)

where J_leak, J_ncx, J_serca, and J_in are defined as in the Monte Carlo approach (see Appendix A ), but

and

are functionals of the probability densities [ρⁱ(c_ds, c_jsr, t)] governed by Eqs. (19), 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) to Eqs. (4)). 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)). 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)), one for each state of the chosen CaRU model (Eq. (3)). Further details of the probability density approach are presented in Appendix B Generalization of the probability density approach,Appendix C Derivation of the univariate probability density approach,Appendix D Numerical scheme for the univariate probability density approach.

Results

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

Figure 2A shows representative Monte Carlo simulations of the minimal whole cell model of EC coupling presented above (Eqs. (1) and Appendix A ). In this simulated voltage-clamp protocol, the holding potential of −80mV is followed by a 20-ms duration test potential to −30, −20, and −10mV (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 −10mV leads to greater Ca²⁺ influx, higher diadic subspace [Ca²⁺], and more Ca²⁺ release than the test potentials of −30 and −20mV. When the test potential is −10mV a 30× “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.

Display large version of this figure

Figure 2

(A) Monte Carlo simulation of the whole cell model exhibits graded release during step depolarization from a holding potential of −80mV to −30, −20, and −10mV (dotted, dot-dashed, and solid lines, respectively). From top to bottom: command voltage, average diadic subspace [Ca²⁺] (

, Eq. (11)), total Ca²⁺ flux via L-type PM Ca²⁺ channels (

, Eqs. (8)), and total Ca²⁺-induced Ca²⁺ release flux (

, Eqs. (8)). The simulation used N=5000 Ca²⁺ release units. (B) Monte Carlo simulations similar to panel A except that the step potential is −10 (solid lines) and +10mV (dotted lines), respectively. Here and below parameters are as in Table 1,Table 2,Table 3.

Figure 2B shows a direct comparison between test potentials of −10 and 10mV. These test potentials result in nearly identical whole cell Ca²⁺ currents (averaged over the duration of the pulse, and 1.4μM/s, respectively). In spite of this, the induced Ca²⁺ release flux is significantly greater when the test potential is −10mV () as opposed to 10mV (21μM/s). This phenomenon occurs because the L-type channel open probability is greater at 10mV than −10mV (Eq. (2)), while the driving force for Ca²⁺ ions is reduced (Eqs. (13)). Although the overall trigger Ca²⁺ flux is nearly the same at these two test potentials, Ca²⁺ release is more effectively induced when the trigger Ca²⁺ is apportioned in larger quantities among a smaller number of diadic subspaces, because the influx that does occur is then more likely to trigger Ca²⁺ sparks. This physiologically realistic aspect of local control during EC coupling is observed in Monte Carlo simulations (see also 19,21), but cannot be reproduced by common pool models 7, nor is it seen in models in which SR Ca²⁺ release depends explicitly on whole-cell Ca²⁺ current (e.g., 16).

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

Display large version of this figure

Figure 3

Solid lines show the dynamics of bulk myoplasmic (c_myo) and network SR (c_nsr) [Ca²⁺] in the whole cell voltage clamp protocol of Fig. 2 with step potential of −10mV (note longer timescale). Dashed lines show the overall myoplasmic (

, Eq. (10)) and network SR (

, Eq. (12)) [Ca²⁺] that include contributions from diadic subspaces and junctional SR, respectively. Note that

is only slightly greater than c_myo and the two traces are not distinguishable.

The dashed line of Fig. 3 shows that the total SR [Ca²⁺] including both network and junctional SR (Eq. (12)) is transiently less than the network SR [Ca²⁺] (), reflecting the fact that for several hundred milliseconds after the voltage pulse junctional SR Ca²⁺ is depleted. While the ratio between the total junctional SR effective volume and the network SR effective volume is , the corresponding ratio between the total diadic subspace volume and the myoplasmic volume is much smaller (). Consequently, the elevated average diadic subspace [Ca²⁺] during the depolarizing voltage step ( as shown in Fig. 2) does not significantly increase the overall myoplasmic [Ca²⁺] ( and the two traces overlap in Fig. 3). On the other hand, depleted junctional SR Ca²⁺ during and after the voltage pulse (, not shown) represents a significant depletion of the overall SR Ca²⁺ content ( in Fig. 3). Although junctional SR depletion develops rapidly after the initiation of the voltage pulse, refilling of these compartments via diffusion of Ca²⁺ from the network SR ( in Eq. (6)) is not complete until ∼400ms after the termination of the voltage pulse (compare solid and dashed lines).

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 −10mV, solid line of Fig. 2). Figure 4A 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 10ms 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).

Display large version of this figure

Figure 4

(A) Dynamics of the diadic subspace (

) and junctional SR (

) Ca²⁺ concentrations associated with a single Ca²⁺ release unit during the voltage clamp protocol of Figure 2 and Figure 3. (B) The dynamics of these local Ca²⁺ concentrations in the (c_ds,c_jsr)-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 [Ca²⁺] given by Eq. (33).

The observation that diadic subspace [Ca²⁺] decreases during the voltage pulse suggests that its dynamics are fast compared to the time evolution of junctional SR [Ca²⁺]. In fact, for the physiologically realistic parameters used in Figure 2 and Figure 3 and Figure 4, the diadic subspace [Ca²⁺] () is well approximated by assuming quasistatic equilibrium with the junctional SR (), bulk myoplasmic (c_myo), and network SR (c_nsr) Ca²⁺ concentrations. Setting the in Eq. (5) and solving for we find that

(33)

where

and

depend on Ca²⁺ release unit state and

and

are functions of plasma membrane voltage defined by

with

as in Eq. (28).

Figure 4B replots the dynamics of the diadic subspace and junctional SR [Ca²⁺] shown in Figure 4A 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 (Figure 2 and Figure 3). Conversely, Fig. 5 presents the state of each of the 5000 CaRUs at a particular moment in time (t=30ms, halfway through the test potential of −10mV). 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).

Display large version of this figure

Figure 5

The open circles are a snapshot at t=30ms of the diadic subspace (

) and junctional SR (

) Ca²⁺ concentrations in the Monte Carlo simulation of Fig. 2. Each of the four central panels corresponds to a particular Ca²⁺ 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) [Ca²⁺] conditioned on CaRU state. Histograms are scaled for clarity and in some cases also truncated (asterisks).

Fig. 5 also shows histograms of the observed distribution of diadic subspace [Ca²⁺] (horizontal) and junctional SR [Ca²⁺] (vertical). The histograms associated with CaRU state clearly indicate that most of these 3387 CaRUs have replete junctional SR (), something that is not obvious from the open circles in the (c_ds, c_jsr)-plane. Similarly, most of the 154 CaRUs in state are associated with replete junctional SR. Conversely, the junctional SR [Ca²⁺] for the 1369 CaRUs in state is broadly distributed with the “average” junctional SR severely depleted (∼100μM). At t=30ms only 90 CaRUs are in state and the distributions of junctional SR [Ca²⁺] and diadic subspace [Ca²⁺] associated with this state are bimodal.

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) 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

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) 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 []. The resulting advection-reaction equations are 26,27,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), but J_efflux* and J_refill* are now functionals of the joint marginal probability densities [],

(44)

(45)

Comparison of probability density and Monte Carlo results

The four histograms presented in Figure 6AD, 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=30ms 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 Figure 6A 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 Figure 6AB, 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 ∼400ms after RyR closure (cf. Fig. 3).

Display large version of this figure

Figure 6

Histograms of the junctional SR Ca²⁺ concentrations (

) at t=30ms in the Monte Carlo simulation of Figure 2 and Figure 3 and Figure 4 and Figure 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 [Ca²⁺] (Eq. (34)) calculated via numerical solution of Eqs. (29). The probability density calculation of the fraction of subunits in each of the four states is denoted by πⁱ (Eq. (46)).

The solid lines of Figure 6AD, show snapshots of the four joint probability densities , , , and as calculated using the probability density approach (t=30ms). These results were obtained by numerically solving Eqs. (29) using the numerical scheme presented in Appendix D (parameters as in Figure 2 and Figure 3 and Figure 4 and Figure 5). Importantly, the entire distribution of junctional SR Ca²⁺ concentrations observed for each CaRU state in the probability density calculation (solid lines) agrees with the corresponding Monte Carlo result (histograms), thereby validating the probability density methodology and our implementation of both approaches. In particular, notice that the fraction of CaRUs in each state given by

(46)

in the probability density calculation is consistent with the Monte Carlo simulation Fig. 5, for example, in Figure 6A

and this corresponds to

in Figure 5A.

While Fig. 6 shows the four marginal probability densities [] for the junctional SR [Ca²⁺] jointly distributed with CaRU state at a particular moment in time, Figure 7A 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 ∼200ms 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–100ms and increases again when t>100ms (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.

Display large version of this figure

Figure 7

Waterfall plot (A) and snapshots (B) of the time evolution of the total probability density for the junctional SR [Ca²⁺] (ρ^T(c_jsr, t) given by Eq.

Article Outline document.write('<a id="close_link" onclick="javascript:document.getElementById(\'fulltext_body\').className=\'outline_closed\';document.getElementById(\'open_link\').style.display=\'block\';return false;" href="#" close="Close Article Outline">X</a>');

Article Information

PubMed

Related Articles

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

Abstract

Introduction

Model formulation

Whole cell model of EC coupling: Monte Carlo formulation

A minimal four-state Ca2+ release unit model

Concentration balance equations

Description of fluxes

Whole cell model of EC coupling: probability density formulation

Summary of model formulation

Results

Representative Monte Carlo simulations

Dynamics of a representative Ca2+ release unit

Dynamics of the population of Ca2+ release units

A univariate probability density formulation for junctional SR [Ca2+]

Comparison of probability density and Monte Carlo results

Article Outline

A minimal four-state Ca²⁺ release unit model

Dynamics of a representative Ca²⁺ release unit

Dynamics of the population of Ca²⁺ release units

A univariate probability density formulation for junctional SR [Ca²⁺]