Article Information

• PDF (253 kb)
• Full Text with Large Figures

PubMed

• Articles by John J. Rice
• Articles by Raimond L. Winslow

Related Articles

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

  Termination of Cardiac Ca Sparks: An Investigative Mathematical Model of Calcium-Induced Calcium Release Biophysical Journal, Volume 83, Issue 1, 1 July 2002, Pages 59-78 Eric A. Sobie, Keith W. Dilly, Jader dos Santos Cruz, W. Jonathan Lederer and M. Saleet Jafri Abstract A Ca spark arises when a cluster of sarcoplasmic reticulum (SR) channels (ryanodine receptors or RyRs) opens to release calcium in a locally regenerative manner. Normally triggered by Ca influx across the sarcolemmal or transverse tubule membrane neighboring the cluster, the Ca spark has been shown to be the elementary Ca signaling event of excitation–contraction coupling in heart muscle. However, the question of how the Ca spark terminates remains a central, unresolved issue. Here we present a new model, “sticky cluster,” of SR Ca release that simulates Ca spark behavior and enables robust Ca spark termination. Two newly documented features of RyR behavior have been incorporated in this otherwise simple model: “coupled gating” and an opening rate that depends on SR lumenal [Ca]. Using a Monte Carlo method, local Ca-induced Ca release from clusters containing between 10 and 100 RyRs is modeled. After release is triggered, Ca flux from RyRs diffuses into the cytosol and binds to intracellular buffers and the fluorescent Ca indicator fluo-3 to produce the model Ca spark. Ca sparks generated by the sticky cluster model resemble those observed experimentally, and Ca spark duration and amplitude are largely insensitive to the number of RyRs in a cluster. As expected from heart cell investigation, the spontaneous Ca spark rate in the model increases with elevated cytosolic or SR lumenal [Ca]. Furthermore, reduction of RyR coupling leads to prolonged model Ca sparks just as treatment with FK506 lengthens Ca sparks in heart cells. This new model of Ca spark behavior provides a “proof of principle” test of a new hypothesis for Ca spark termination and reproduces critical features of Ca sparks observed experimentally. Abstract | Full Text | PDF (936 kb)

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

  Moment Closure for Local Control Models of Calcium-Induced Calcium Release in Cardiac Myocytes Biophysical Journal, Volume 95, Issue 4, 15 August 2008, Pages 1689-1703 George S.B. Williams, Marco A. Huertas, Eric A. Sobie, M. Saleet Jafri and Gregory D. Smith Abstract In prior work, we introduced a probability density approach to modeling local control of Ca-induced Ca release in cardiac myocytes, where we derived coupled advection-reaction equations for the time-dependent bivariate probability density of subsarcolemmal subspace and junctional sarcoplasmic reticulum (SR) [Ca] conditioned on Ca release unit (CaRU) state. When coupled to ordinary differential equations (ODEs) for the bulk myoplasmic and network SR [Ca], a realistic but minimal model of cardiac excitation-contraction coupling was produced that avoids the computationally demanding task of resolving spatial aspects of global Ca signaling, while accurately representing heterogeneous local Ca signals in a population of diadic subspaces and junctional SR depletion domains. Here we introduce a computationally efficient method for simulating such whole cell models when the dynamics of subspace [Ca] are much faster than those of junctional SR [Ca]. The method begins with the derivation of a system of ODEs describing the time-evolution of the moments of the univariate probability density functions for junctional SR [Ca] jointly distributed with CaRU state. This open system of ODEs is then closed using an algebraic relationship that expresses the third moment of junctional SR [Ca] in terms of the first and second moments. In simulated voltage-clamp protocols using 12-state CaRUs that respond to the dynamics of both subspace and junctional SR [Ca], this moment-closure approach to simulating local control of excitation-contraction coupling produces high-gain Ca release that is graded with changes in membrane potential, a phenomenon not exhibited by common pool models. Benchmark simulations indicate that the moment-closure approach is nearly 10,000-times more computationally efficient than corresponding Monte Carlo simulations while leading to nearly identical results. We conclude by applying the moment-closure approach to study the restitution of Ca-induced Ca release during simulated two-pulse voltage-clamp protocols. Abstract | Full Text | PDF (568 kb)

• A Probability Density Approach to Modeling Local Control of ...

  A Probability Density Approach to Modeling Local Control of Calcium-Induced Calcium Release in Cardiac Myocytes Biophysical Journal, Volume 92, Issue 7, 1 April 2007, Pages 2311-2328 George S.B. Williams, Marco A. Huertas, Eric A. Sobie, M. Saleet Jafri and Gregory D. Smith 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. Abstract | Full Text | PDF (706 kb)

• …more

Copyright © 1999 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 77, Issue 4, 1871-1884, 1 October 1999

doi:10.1016/S0006-3495(99)77030-X

Biophysical Theory and Modeling

Modeling Gain and Gradedness of Ca²⁺ Release in the Functional Unit of the Cardiac Diadic Space

John J. Rice, M. Saleet Jafri^,   and Raimond L. Winslow

Department of Biomedical Engineering and Center for Computational Medicine and Biology, The Johns Hopkins University School of Medicine, Baltimore, Maryland 21205 USA

Address reprint requests to Dr. M. Saleet Jafri, Traylor Research Building, Room 412, 720 Rutland Avenue, Baltimore, MD 21205. Tel.: 410-502-5091; Fax: 410-614-0166.

Calcium-induced calcium release (CICR) in cardiac muscle exhibits both gradedness and high gain. Gradedness refers to the observation that sarcoplasmic reticulum (SR) Ca²⁺ release is proportional to the influx of trigger Ca²⁺ (Beuckelmann and Wier, 1988), whereas high gain refers to the observation that Ca²⁺ release from SR is significantly larger than the trigger influx (Fabiato, 1985a). The paradox, as described by Stern, 1992 is that the positive feedback inherent in such high-gain systems produces all-or-none rather than graded Ca²⁺ release. Such behavior is predicted for all models in which sarcolemmal Ca²⁺ influx and SR Ca²⁺ release are directed into a single compartment (referred to as “common pool” models).

Our previous model of cardiac calcium handling incorporated an improved description of the L-type Ca²⁺ channels and CICR release from ryanodine receptor (RyR) channels exhibiting adaptation. Because these channels, as well as the diadic space between the channels, are each represented as ensemble averages, our previous work is an example of a common pool model (Jafri et al). While this model reproduces important frequency-dependent aspects of cardiac Ca²⁺ cycling and high gain, it also exhibits all-or-none rather than graded Ca²⁺ release. We hypothesize that much different behavior will result from a stochastic implementation in which individual functional release units (FRUs), each consisting of dihydropyridine (DHPR) and RyR channels interacting via a local functional unit subspace. This hypothesis assumes that a variable number of independent FRUs can be recruited via a “local control” mechanism to produce graded Ca²⁺ release (Stern, 1992,Cannell et al,Lopez-Lopez et al). This assumption is consistent with the spatially and temporally localized SR Ca²⁺ release events, known as Ca²⁺ sparks, that are thought to be the unitary CICR events in cardiac cells (Lipp and Niggli, 1994,Cheng et al). We assume that single voltage-activated DHPR channels can open nearby RyRs via localized increases of [Ca²⁺] (Cannell et al,Lopez-Lopez et al). This system is of sufficiently low order to be run with stochastic (Monte Carlo) simulations, and allows for systematic parameter variation. The model shows robust graded Ca²⁺ release and produces behavior consistent with that reported for Ca²⁺ sparks.

Fig. 1 shows the model of the functional release unit. The model consists of one DHPR and eight RyR channels communicating via a functional unit subspace. Experimental studies have suggested that a single DHPR opening activates a single functional release unit (Sham et al,Cannell et al). The 8:1 RyR/DHPR stoichiometry used is similar to the experimentally determined value of 7.3 (Bers and Stiffel, 1993). Calcium enters the subspace by two pathways: across the sarcolemma via a single DHPR and from the SR via any of eight RyR channels. These two Ca²⁺ fluxes are labeled J_DHPR and J_RyR, respectively. The subspace is considered to be a single compartment with no spatial Ca²⁺ gradients. Therefore, each RyR channel will sense the same Ca²⁺ concentration ([Ca²⁺]_SS) within the subspace. This simplifying assumption is based on previous three-dimensional models of Ca²⁺ diffusion near the channel pore. They suggest that [Ca²⁺] a short distance away from the channel parallel to the membrane (tens of microseconds) almost instantaneously reaches a uniform spatial profile (Langer and Peskoff, 1996,Soeller and Cannell, 1997,Keizer and Smith, 1998). For simplicity, the myoplasmic Ca²⁺ concentration ([Ca²⁺]_myo) is assumed to be a fixed constant (0.1μM).

Display large version of this figure

Figure 1

Schematic diagram of the functional release unit model. Ca²⁺ influx through the DHPR (J_DHPR) enters the functional unit subspace (dashed box). This rise in local [Ca²⁺] fosters the opening of one or more of the eight RyR, producing additional influx (J_RyR) into the subspace. In this model, the subspace is considered to be a single compartment, hence there is no spatial gradient across this space. Efflux from the subspace (J_xfer) occurs via simple diffusion, as does the flux to refill junctional SR (J_tr).

We also assume that each FRU behaves independently of others based on experimental findings (Sham et al,Cannell et al,Cheng et al). Cheng and co-workers observed that Ca²⁺ sparks become propagating Ca²⁺ waves only under Ca²⁺ overload conditions, but not under normal conditions, suggesting that activation does not typically spread from one functional unit to the next. Sparks are also thought to consist of Ca²⁺ release from between 6 and 20 RyRs, a finding that is consistent with the RyR/DHPR ratio used in this study. Typically, during a contraction the activation of each functional unit is spread by Ca²⁺ entry via the DHPR during depolarization of the sarcolemma. Ca²⁺ flux from the subspace to myoplasm (J_xfer) is computed as the [Ca²⁺] gradient divided by the relaxation time constant (τ_xfer):

(1)

(2)

The DHPR channel is represented by a mode-switching model that was developed previously (Jafri et al). A state diagram of the mode switching and voltage activation is shown in Figure 2A. Transition rates are given in Table 1. The upper and lower rows of states comprise Mode Normal and Mode Ca, respectively. The channel is assumed to be composed of four independent subunits that can each close the channel. This dictates five closed states (C₀-C₄) on the top row and a mirror set of closed states (C_Ca0-C_Ca4) on the bottom row. The proportionality of the forward and reverse rates between the closed states is dictated by the four-way symmetry assumed for channel subunits. Voltage-dependent activation is incorporated through the rate constants α and β, which are increasing and decreasing functions of voltage, respectively. When in the rightmost closed states (C₄ or C_Ca4), there are voltage-independent transitions to the open states (O or O_Ca). Note that f′ is 500 times slower than f, so that openings are rare in Mode Ca, effectively inactivating the channel.

Display large version of this figure

Figure 2

State diagrams for the DHPR and RyR channel Markov models. (A) A state diagram of the mode switching and voltage activation of DHPR channel model. The upper row of states comprises Mode Normal and the lower row comprises Mode Ca. The channel is composed for four independent subunits, each of which can close the channel. The corresponding states are C_n, where n is the number of permissive subunits. With depolarization, the channels undergo transitions from left to right. With four permissive subunits (C₄), there is a voltage-independent transition to the conducting state O. With elevation of [Ca²⁺] the channel transition occurs from the Mode Normal (top) to Mode Ca (bottom). f ≪ f′ so that transitions into O_Ca are rare, which effectively inactivates the channel in Mode Ca. (B) A state model for the RyR with adaptation. The states P_C1 and P_C2 are the resting states, P_O3 and P_O4 are the open states, and P_C5 and P_C6 are the adapted or refractory states. The transitions governed by the rates k_1,2, k_2,3, k_3,4, k_5,4, k_5,6, and k_2,5 are all Ca²⁺-dependent (see Table 2 for values of the rate constants).

Table 1 DHPR transition probabilities

	γ=0.1875 [CaSS] s−1	a=2.0s−1
	α=0.4 e(V+6)/25 s−1	b=2.0s−1
	β=0.05 e−(V+6)/29 s−1	f=300.0s−1
	α′=αa s−1	g=2000.0s−1
	β′=(β/b)s−1	f′=5.0s−1
	ω=10.0s−1	g′=7000.0s−1
	y∞=0.02+
	τy=0.02+

Transitions to Mode Ca are controlled by γ, which is a function of Ca²⁺. Moving rightward in Figure 2A, there are incremental increases in the multiplier of γ and the divisor on ω. The effect of this is to greatly increase the transition rate to Mode Ca at high voltages when the channel is opening. The close symmetry between Mode Normal and Mode Ca closed states and similarity of rates (i.e., g versus g′) is dictated by the experimental finding that gating currents are very similar in the Mode Normal and Mode Ca cases, and by thermodynamic microscopic reversibility. Microscopic reversibility requires that for each cycle the product of rates is equal, whether taken in the clockwise or the counterclockwise direction. Voltage-dependent inactivation is modeled as a Hodgkin-Huxley-type gate (y). This gate can inactivate the channel independently of the states discussed above. The equations modeling the DHPR are provided in Table 2.

Ca²⁺ flux through the DHPR channel (J_DHPR) is computed as

(3)

where Ī_DHPR is open channel current (μA μF⁻¹) given by

(4)

DHPR_open is 1 for states O or O_Ca, and is 0 otherwise; F is Faraday's constant, and V_SS is the volume of the FRU subspace. The number 0.341 is the activity coefficient of Ca²⁺ at the outer opening of the cell. The number 0.001 reflects the product of the activity coefficient of Ca²⁺ at the mouth of the channel in the FRU subspace and [Ca²⁺]_SS.

The RyR channel is represented using a model developed by Keizer and Smith, 1998. This model was developed to replicate open and dwell times of isolated RyR channels in vitro and in vivo, as well as measured peak and plateau open probabilities with Ca²⁺ or cesium (Cs²⁺) as the charge carrier. The latter measurements describe the adaptive behavior of the RyR channel. As originally described, adaptation is a property of the RyR where, after rapid activation by a step increase in [Ca²⁺], the channel undergoes a slow spontaneous decrease in open probability (Györke and Fill, 1993). Closing of the RyR has also been attributed to inactivation (Fabiato, 1985b,Sham et al). In isolated bilayers, adaptation occurs within milliseconds, while inactivation occurs within a few seconds (Keizer and Levine, 1996). Hence, in this model we are interested in the effects of adaptation on RyR Ca²⁺ release.

A state diagram of the RyR model is shown in Figure 2B. The state P₀₃ and P_O4 are open states whereas states P_C1, P_C2, P_C5, and P_C6 are closed states. Briefly, P_C1 and P_C2 represent the initial closed states. With an increase in [Ca²⁺], the channel opens, moving to states P₀₃ and P_O4. The other states, P_C5, and P_C6, represent the adapted or refractory state. Transition rates between the states are given by k_x,y, and are provided in Table 3. In the original model (Keizer and Smith, 1998), the charge carrier could be Ca²⁺ or Cs²⁺, so that some transition rates depended on [Ca²⁺] in the microdomain around the channel while others depended on the bulk myoplasmic [Ca²⁺]. We assume the charge carrier is Ca²⁺ so all rates depend on the [Ca²⁺] in the microdomain. Furthermore, we have modified the rates so that 1) they are scaled to depend on FRU subspace Ca²⁺ and not on bulk myoplasmic Ca²⁺, and 2) they are saturating functions of Ca²⁺ following Michaelis-Menten kinetics. The new description was derived to match the original description over low levels of [Ca²⁺], but prevents extremely large, and likely unrealistic, transition rates at high [Ca²⁺]. Such high rates would require extremely small time steps that greatly increase computation time. The results of both these modifications are shown in Table 2.

Table 3 Parameters

	VSS=2.03×10−12μL	[BSL] = 1124.0μM
	VJSR=1.05×10−10μL	KBSL=8.7μM
	[Ca2+]myo=0.1μM	[CSQNtotal] = 15.0mM
	[Ca2+]NSR=800μM	KCSQN=0.8mM
	[Ca2+]o=1.8mM	F=96500 coul (mol e−)−1
	τtr=0.005ms	T=310K
	τxfer=0.0007ms	R=8.314J mol−1K−1
	[BSR]=47.0μM	= 3960.0s−1
	KBSR=0.87μM	=1.51×10−11cmms−1

The Ca²⁺ flux through the RyR channels (J_RyR) is computed as

(5)

where is the channel permeability to Ca²⁺, RyR_openⁱ is 1 when the ith channel is in state P₀₃ or P_O4 and 0 otherwise, and ([Ca²⁺]_SS−[Ca²⁺]_JSR) is the driving force from the JSR to the subspace.

The full Ca²⁺ balance equations for the subspace and the JSR are

(6)

and

(7)

where

(8)

and

(9)

describe the buffering of Ca²⁺ using the rapid buffer approximation (

Monte Carlo simulations are run using standard methods for Markov processes (Keizer, 1987). Unless noted otherwise, all simulations use a voltage clamp protocol. Specifically, the model is allowed a period of time (0.1s) to reach steady state at a holding potential of −80mV, followed by a voltage step to a test potential for 0.2s, and then returned to the −80mV holding potential. Each trial is repeated 500 times before changing the voltage step (typically increased from −50 and +50mV in 5-mV increments). Results are averaged over the 500 independent trials at each voltage, hence could be considered the ensemble behavior of 500 independent functional units. This number is probably less than the total number of functional units in a cell, which is estimated to be in the range of 10³-10⁴ (Isenberg, 1995). This number was selected to yield a reasonable compromise between minimization of run time and reduction of variability in the results.

Sample Monte Carlo results are shown in Fig. 3 for a single trial with a voltage step to 0mV. The peak fluxes are shown for the DHPR channel (solid line) and the sum of the eight RyRs (dashed line) in Figure 3A. After a voltage transition to 0mV at time 0.1s, the DHPR might open (Figure 3A). In the simulation shown, the DHPR opens only once. Often there will be no openings and sometimes more than one opening. The DHPR flux displays a consistent level of peak flux for a given clamp voltage, because we assume a single channel with constant flux that depends only the clamp voltage. In contrast, the RyR show a variable level of flux because 1) a subset of eight channels is open at any given time; and 2) the driving force for Ca²⁺ flux depends on [Ca²⁺]_JSR, which varies with time. Note that RyR flux is largest at the beginning of the voltage step, reflecting a relatively large driving force and higher open probability early in the Ca²⁺ release event. This effect can also be seen in Figure 3B, where the flux is integrated over time to produce the total integrated flux of Ca²⁺ into the subspace from each of the two possible sources. Figure 3CD show the corresponding changes in Ca²⁺ concentration in the subspace ([Ca²⁺]_SS) and the JSR ([Ca²⁺]_JSR), respectively. The peak [Ca²⁺]_SS reaches a value of ∼22μM. The local JSR depletes from 800 to ∼400μM.

Display large version of this figure

Figure 3

Sample results for a single functional unit driven by a voltage step to 0mV from 0.1 to 0.3s. (A) Peak fluxes are shown for the DHPR channel (solid line) and the sum of the eight RyRs (dashed line). (B) Integrated fluxes are shown for the DHPR channel (solid line) and the sum of the eight RyRs (dashed line). (C) The changes in Ca²⁺ concentration in the subspace ([Ca²⁺]_SS). (D) The changes in JSR Ca²⁺ concentration ([Ca²⁺]_JSR) indicates the degree to which SR empties during a Ca²⁺ release event.

While Fig. 3 shows results for a single trial, Fig. 4 shows similar data averaged over 500 independent runs. The average peak and integrated flux through DHPR (lower line) and RyRs (upper line) are shown in Figure 4AB, respectively. In both cases, the magnitude of the averaged data is somewhat lower than the corresponding single channel data because of the occurrence of runs in which the DHPR channel does not open, thus producing no RyR openings. Figure 4CD show the corresponding average changes in [Ca²⁺]_SS and [Ca²⁺]_JSR. The change in [Ca²⁺]_JSR from a resting value of 800 to 562μM was a decrease of 38%. A similar decrease is observed experimentally as cardiac myocytes retain 35% of their resting SR Ca²⁺ during a contraction under normal SR loading conditions (Bassani et al).

Display large version of this figure

Figure 4

Average results for 500 single runs (same protocol as used in Fig. 3 for a single functional unit driven by a voltage step to 0mV from 0.1 to 0.3s. (A) Average peak fluxes are shown for the DHPR channel and the sum of the eight RyRs. (B) Average integrated fluxes are shown for the DHPR channel (solid line) and the sum of the eight RyRs (dashed line). (C) The changes in average Ca²⁺ concentration in the subspace ([Ca²⁺]_SS). (D) The changes in average JSR Ca²⁺ concentration ([Ca²⁺]_JSR) indicate the degree to which SR empties during a Ca²⁺ release event.

Again, the data are smoother and show smaller changes in magnitude when compared to the corresponding data for a single functional unit in Fig. 3. These results demonstrate that the model produces positive gain, i.e., average integrated flux through the RyRs is ∼10 times larger than that through the DHPR.

The next question to address is whether Ca²⁺ release is graded with influx of trigger Ca²⁺ through the DHPR. To test gradedness, the step potential is varied from −50 to +50mV in 5-mV increments. When plotted against clamp voltage, the averaged DHPR fluxes produce bell-shaped profiles for both integrated (Figure 5A, solid line) and peak fluxes (Figure 5C, solid line). The peak of the DHPR flux is ∼−10mV, similar to that measured experimentally. Likewise, RyR flux (dashed lines) shows similar bell-shaped profiles, but are much larger in amplitude. These data clearly show that the model can reproduce both high gain and graded Ca²⁺ release that is observed in cardiac myocytes. The gain or amplification factor is defined as the ratio of the RyR flux to the DHPR flux. It shows the amount of amplification of the trigger Ca²⁺ by CICR. The gain for the integrated fluxes is shown in Figure 5B and the gain for the peak flux is shown in Figure 5D. While the gain shown is in the range observed physiologically (Wier et al), the shape differs from that seen in the experiments. Possible explanations for this will be described in the Discussion.

Display large version of this figure

Figure 5

Model exhibits positive gain and graded Ca²⁺ release. To modulate the amount of trigger Ca²⁺, voltage pulses are changed in 5-mV steps from −50 to 50mV. The data are averaged from 500 runs. (A) Integrated flux through the DHPR (solid line) shows a bell-shaped profile, similar to those measured experimentally. The integrated RyR Ca²⁺ release (dashed line) also has a bell-shaped profile, which shows gradedness and is much larger, that demonstrates positive gain. (B) The gain or amplification factor for the integrated fluxes. (C) Peak flux shows very similar behavior, although there are slight differences in the shapes of the curves when compared to the integrated flux data. (D) The gain for the peak fluxes.

In DHPR channels, separate features produce the rising and falling phases of the bell-shaped flux profile. The rising phase is produced by the increasing open probability and results from voltage-sensitive channel activation. The declining phase reflects changes in open-channel current. The two phases of the RyR flux curves also reflect this dual modulation. This is demonstrated in the next set of simulations. In Figure 6A, the activation function of the DHPR channel is shifted by −10 and +10mV, producing shifts of the rising edge of the peak DHPR flux functions (seen by the three low-amplitude curves in Figure 6A). The rising phases of the RyR peak flux functions show corresponding shifts (seen by the three higher amplitude curves in Figure 6A). Note that there are only small changes in the declining phases because the open probability of the DHPR is close to its saturating level. The declining phase predominantly reflects the open-channel I-V relation of the DHPR functions. As shown in Figure 6B, shifting the open-channel I-V relation by −10 and +10mV produces corresponding changes in both DHPR and RyR flux functions. These results indicate that properties of the RyR release flux are closely coupled to processes governing activation and permeation of the DHPR.

Display large version of this figure

Figure 6

The rising phase of the bell-shaped profiles reflects voltage-dependent activation of the DHPR channel, whereas the falling phase reflects the open channel I-V relations. To modulate the amount of trigger Ca²⁺, voltage pulses are changed in 5-mV steps from −50 to 50mV. The data are averaged from 500 runs. (A) The activation function of the DHPR channel is shifted by −10 (dashed line) and +10mV (dotted line). The rising phase of the RyR average peak flux curves shows corresponding shifts, showing that RyR Ca²⁺ release closely corresponds to DHPR open probability on the rising phase. (B) The I-V relations of the DHPR are shifted by −10 (dashed line) and +10mV (dotted line). The falling phase of the RyR average peak flux curves show corresponding shifts, showing that RyR Ca²⁺ release closely corresponds to DHPR I-V relations on the falling phase.

The next set of simulations addresses the role of adaptation of RyR Ca²⁺ release. RyR adaptation rate is first increased and then decreased by a factor of 5 by adjusting the transition rates k_4,5 and k_2,5 (see Figure 2B). A slower adaptation rate increases peak RyR flux and [Ca²⁺] in the subspace (Figure 7AB, dotted lines). Conversely, the effect of more rapid adaptation is to decrease both peak flux and [Ca²⁺]_SS (Figure 7AB, dashed lines). It is also clear that rapid adaptation produces a model with a significantly decreased gain that is smaller than that observed experimentally (Wier et al). This finding suggests that measurements of gain may help to limit estimates of RyR adaptation rates.

Display large version of this figure

Figure 7

The effects of modifying adaptation rate on RyR Ca²⁺ release. The transition rates into the adapted state (k_2,5 and k_4,5 in Figure 2B) are increased or decreased by a factor of 10 to increase or decrease adaptation rate, respectively. Data are gathered from voltage pulses from −50 to 50mV in 5-mV increments, and are averaged from 500 runs at each voltage. (A) Slow adaptation rate (dotted line) produces a larger peak flux than control conditions (solid line), the same as Figure 5C). Increased adaptation rate decreases the peak flux. (B) The average [Ca²⁺]_SS transients at 0mV show that adaptation rate effects both the magnitude and temporal profile of RyR Ca²⁺ release. (C) The average [Ca²⁺]_JSR lines show that the higher rate of adaptation (dashed line) decreases the degree of SR unloading during Ca²⁺ release, whereas decreased rate (dotted line) produces greater SR unloading.

The effect of RyR adaptation on JSR Ca²⁺ concentration is examined in Figure 7C. Slowed adaptation produces a larger SR Ca²⁺ release and lower [Ca²⁺]_JSR during the course of the Ca²⁺ release (dotted line). A more rapid adaptation rate produces a smaller SR Ca²⁺ release and higher [Ca²⁺]_JSR during the course of the Ca²⁺ release (dashed line). Therefore, in this model, RyR adaptation rate has an effect on the level of SR emptying and contributes to termination of SR Ca²⁺ release.

Considerable experimental evidence suggests that SR Ca²⁺ release is positively correlated with the Ca²⁺ load (Bers, 1991,Bassani et al; Janczewski et al). For example, by decreasing the SR Ca²⁺ content by 56%, its control value reduces the RyR Ca²⁺ release by 52% (Janczewski et al). The effect of altered SR Ca²⁺ load is tested in the model by changing the initial Ca²⁺ concentration of JSR ([Ca²⁺]_JSR) and the concentration of the store that refills SR ([Ca²⁺]_NSR). When SR Ca²⁺ load is doubled from the control value (1600μM vs. 800μM), the peak flux increases over the entire range of voltages tested (Figure 8A, dashed line). Similarly, when SR Ca²⁺ load is halved from the control value (400μM vs. 800μM), the peak flux decreases (Figure 8A, dashed line).

Display large version of this figure

Figure 8

The effects of modifying SR Ca²⁺ load on RyR Ca²⁺ release. The [Ca²⁺] of JSR and the [Ca²⁺] of the refilling store are doubled (1600μM) (dashed line) or halved (400μM) (dotted line) to effectively increase or decrease the SR load. Data are from voltage pulses from −50 to 50mV in 5-mV increments and are averaged from 500 runs at each voltage. (A) Increased SR load (dashed line) produces a larger peak flux than control conditions (solid line), the same as Figure 5C). (B) The average [Ca²⁺]_SS transients at 0mV show that SR load effects both the magnitude and temporal profile of RyR Ca²⁺ release. Doubling the SR load (dashed line) produces elevated sustained Ca²⁺ release during the voltage pulse (0.1–0.3s) and regenerative RyR Ca²⁺ release after the pulse when DHPR influx has ceased. (C) The average [Ca²⁺]_JSR lines show the degree of SR unloading during and after the voltage pulse.

The change in SR load also produces large changes in the temporal behavior of SR Ca²⁺ release. As shown in Figure 8B, reduced SR load (dotted line) produces a smaller, shorter peak in [Ca²⁺]_SS with a smaller sustained component after the Ca²⁺ release. In contrast, increased SR load produces a larger magnitude of Ca²⁺ release with an elevated sustained component. Ca²⁺ release also continues after the voltage step (see dashed line after the voltage step from 0.1–0.3s in Figure 8B). This phenomenon occurs only when SR Ca²⁺ load is elevated so that release become regenerative as RyRs continue to activate each other (i.e., the cluster bomb effect as described by Stern, 1992).

In contrast to regenerative RyR openings, the data of Fig. 8 demonstrate that model RyR Ca²⁺ release closely follows the DHPR openings. It has been pointed out elsewhere (Callewaert, 1992,Stern, 1992) that such a system produces high gain and graded Ca²⁺ release. However, the data of Figure 8B also show that under at least some conditions (i.e., high SR load), RyR Ca²⁺ release can be sustained in the absence of influx of trigger Ca²⁺.

In the next set of simulations, we seek to determine the functional coupling between the DHPR and RyRs as compared to coupling between adjacent RyRs. If the former is much stronger, then RyRs may simply follow the trigger influx. If the latter is strong, then the system is more likely to tend toward self-regenerating RyR Ca²⁺ release. To address this issue, the stochastic model is modified so that either a single DHPR or single RyR channel is held open to drive the functional unit in the absence of any other trigger influx of Ca²⁺. An important variable in such simulations is the open duration. Therefore, an initial step is to determine an appropriate range of channel open times for use in the model.

Fig. 9 shows the channel open times plotted as a function of initial time of the opening event. The model was run using the same standard conditions as in Figure 3 and Figure 4 and Figure 5. In Figure 9A, the DHPR opens times are plotted during the 0.1–0.3-ms voltage step to 0mV. Most points are below 3ms, and the distribution of points shows that although there are more frequent openings early in the voltage step, there is no apparent dependence of open time on the time of opening. To see that this is the case, the data points within 1-ms time intervals were averaged, and the result is plotted as a gray line in Figure 9A. The line has slope near zero, demonstrating that the DHPR mean open time is constant at ∼0.5ms. This corresponds well with the value of 0.45ms measured experimentally by Rose and co-workers, 1992. In Figure 9B, the RyR open times show some time dependence with a peak mean open time of ∼6.7ms early during the voltage pulse and mean open time of ∼2.8ms later in the pulse. These values are in agreement with the ranges of values from 1.5 to 7.1ms measured experimentally for control conditions (Lukyanenko et al,Eager and Dulhunty, 1998). For 500 runs, the 14,520 RyR openings (Figure 9B) are much more numerous than 1567 DHPR openings (Figure 9A).

Display large version of this figure

Figure 9

The DHPR and RyR open times during a voltage pulse to 0mV under default conditions. The channel open times (ordinate) plotted as a function of initial time of the opening event (abscissa). (A) DHPR open times show little dependence on initial time of opening, as shown by the average for 1-ms bins (white-black line). The mean open time is 0.50±0.51ms for 1567 points collected from 500 runs. (B) RyR open times show some time-dependence with a peak average open time of ∼6.7ms early during the voltage pulse, and average around 2.8ms (white line) later in the clamp. The mean open time is 2.8±3.8ms for 14,520 points collected from 500 runs.

Fig. 10 shows model results when a single DHPR or RyR channel is held open. The average Ca²⁺ transients for 500 runs are shown in Figure 10A for a DHPR open time of 0.5ms (black line) and 16.0ms (gray line). The amplitudes are similar with a slight variation of the time course. This compares well with the experimental results of Sham et al in which the mean open time of the DHPR is increased from 0.27ms to 15.9ms by the addition of the agonist FPL64176. The ensuing Ca²⁺ transients caused by Ca²⁺ release from the SR are similar in amplitude, with the Ca²⁺ transient with the agonist being slightly longer in duration. When a single RyR is held open for 0.5ms (black line) and 16.0ms (gray line), there is virtually no difference between the FRU subspace Ca²⁺ transients. These simulations are run with the default initial SR Ca²⁺ load ([Ca²⁺]_JSR=[Ca²⁺]_NSR=800μM) so the SR load is similar to those during initial portions of the voltage step in previous simulations (Figure 3 and Figure 4 and Figure 5). The degree of SR unloading is also similar to previous simulations. For the 5- and 16-ms single channel openings, the [Ca²⁺]_JSR falls to 370 and 394μM, respectively, for the DHPR case and 376 and 394μM for RyR case (data not shown). [Ca²⁺]_JSR can be clamped at 800μM so that no SR emptying occurs. The corresponding Ca²⁺ transients are of similar magnitude, but are over 10 times longer (>0.6s, data not shown). This finding suggests that SR unloading plays an important role in the termination of SR Ca²⁺ release.

Display large version of this figure

Figure 10

The RyR Ca²⁺ release elicited from a single DHPR or RyR opening. (A) The average Ca²⁺ transients for 200 runs with default initial SR Ca²⁺ load ([Ca²⁺]_JSR=800μM). The response for a 0.5ms (black line) and a 16ms (gray line) DHPR opening are similar, although the latter is slightly larger and has a small shoulder as a result of the DHPR Ca²⁺ influx. (B) The Ca²⁺ transients for a single a 0.5ms (black line) or 16ms (gray line) RyR opening are similar, although slightly smaller that those produced by a single DHPR opening. (C) RyR Ca²⁺ release is independent of the open duration of a single DHPR, but highly dependent on SR load. The [Ca²⁺] of JSR and the [Ca²⁺] of the refilling store are set to 800, 600, 400, and 200μM (as labeled) to vary the SR load. (D) RyR Ca²⁺ release shows some dependence on the open duration of a single RyR for short open times. As in (C) above, RyR Ca²⁺ release is very dependent on SR load.

The above data show that both single DHPR and single RyR openings can trigger Ca²⁺ release from adjacent RyRs. This suggests that termination of SR Ca²⁺ release may be quite difficult. However, as seen previously in Figure 4A, RyR Ca²⁺ release stops shortly after DHPR influx ceases. A possible explanation is that RyR open times tend to become shorter during later phases of the voltage pulse (Figure 9B), which could produce less activation of neighboring RyRs. However, the possibility that RyR open time alone can lead to less regenerative release is unlikely, given that a 0.5-ms RyR opening still produces a Ca²⁺ transient similar to a 16.0-ms RyR opening (Figure 10B). A short RyR opening can still be very effective in triggering release from its neighbors in the model.

A more important feature for termination of SR Ca²⁺ release is the depletion of SR Ca²⁺ that occurs during the voltage step. This is demonstrated by repeating the single channel opening experiments for a range of open times and for a range of SR loads. Results are shown in Figure 10CD. The peak RyR Ca²⁺ flux for different values of initial [Ca]_JSR are shown. In Figure 10C, the DHPR was held open for a range of open times from 0.1ms to 50.0ms. The peak RyR Ca²⁺ flux is invariant of DHPR open time, as observed by Sham et al, but depends on the initial SR Ca²⁺ concentration. The results show that integrated RyR Ca²⁺ flux depends mainly on the [Ca]_JSR, with little dependence on open duration. In contrast, DHPR open duration has almost no effect on integrated RyR flux, except at very short open durations. The reason is that as the duration of the DHPR opening increases, the RyR tend to adapt so that further RyR Ca²⁺ release is not produced. When the adaptation rate is increased by a factor of 10, integrated RyR flux decreases to ∼10% of the control value, but again with little dependence on DHPR open duration (data not shown). Decreasing adaptation rate by a factor of 10 increases integrated RyR flux by 50%. Again, integrated Ca²⁺ release flux shows little dependence on DHPR open duration (data not shown).

In contrast to the DHPR data, RyR open time can strongly influence peak RyR flux. As shown in Figure 10D, integrated RyR flux strongly depends on the [Ca]_JSR. The higher [Ca]_JSR, the greater the peak flux. Integrated RyR flux also increases monotonically with RyR open duration.

The final set of simulations studies the effects of altered FRU subspace geometry on graded Ca²⁺ release. These simulations are motivated by experiments by Gomez and co-workers showing SR Ca²⁺ release decreases during congestive heart failure (CHF; Gomez et al). The researchers observed that the ability of the DHPR to activate SR Ca²⁺ release was reduced in the hearts of CHF rats. The suggested that a possible mechanism for this is that the DHPRs and RyRs become functionally uncoupled, possibly as a result of altered geometry in diadic space. Here, we seek to uncouple DHPR and RyRs in the model by one of two ways, increasing FRU subspace volume or increasing the efflux of Ca²⁺ out of the FRU subspace.

Fig. 11 shows the effects of altered functional unit subspace volume. Doubling the FRU subspace volume (Figure 11A, dashed line) produces a larger peak flux than control conditions (solid line, same as Figure 5C), with the peak shifted to lower step potentials. Hence, the larger subspace volume can produce larger release flux at low voltages, but decreases the effectiveness of the DHPR trigger at high voltages. Decreasing the FRU subspace volume by half increases the peak flux and shifts it to higher step potentials (Figure 11A, dotted line) compared to the control (solid line). In this case, the small DHPR trigger Ca²⁺ influx seen at high voltages is more effective at triggering RyR Ca²⁺ release than in the control. However, a smaller subspace produces higher subspace [Ca²⁺] reducing release at lower voltages (Figure 11B, dotted line), and the depletion of Ca²⁺ from the SR is less for the smaller subspace volume (not shown).

Display large version of this figure

Figure 11

The effects of modifying subspace volume on RyR Ca²⁺ release. The volume of the subspace is doubled (dashed line) or halved (dotted line). Data are from voltage pulses from −50 to 50mV in 5-mV increments and is averaged from 500 runs at each voltage. (A) Larger subspace (dashed line) produces a larger peak flux than control conditions (solid line), the same as Figure 5C), but only at lower voltages. At higher voltages, the trigger influx through the DHPR channel is insufficient to trigger Ca²⁺ release. The smaller subspace volume (dotted line) produces greater Ca²⁺ release at high voltage because the small amount of trigger is more effective. However, overall gain is lower because RyR adaptation is increased. (B) The average [Ca²⁺]_SS transients at 0mV show that subspace volume effects both the magnitude and temporal profile of RyR Ca²⁺ release. Halving the subspace (dotted line) produces elevated sustained Ca²⁺ release during the voltage pulse (0.1 to 0.3s) and a small amount regenerative RyR Ca²⁺ release after the pulse when DHPR influx has ceased. Note that doubling the subspace increased peak flux at 0mV (A), but produced a smaller change in [Ca²⁺]_SS given the larger volume. (C) Doubling the transfer rate out of the subspace (τ_xfer) results decreases RyR Ca²⁺ release at high depolarizations (dotted line). Decreasing τ_xfer by half increases the gain of RyR Ca²⁺ at higher voltages. (D) The average [Ca²⁺]_SS transients at 0mV show a higher peak value when τ_xfer is decreased (dashed line) and a lower peak value when τ_xfer is increased (dotted line).

The second method involves modification of the transfer rate out of the FRU subspace by altering the time constant τ_xfer (transfer rate is inversely proportional to J_xfer). When the transfer rate is doubled, the right side of the graded release curve (which is determined by the I-V relation) is shifted to the left (Figure 11C, dotted line). Similar to the increased subspace volume case, the doubled transfer rate decreases the effectiveness of the DHPR trigger at high voltages (Figure 11D, dotted line). Reducing the transfer rate by half causes the right side of the graded release curve to be shifted to the right (Figure 11C, dashed line). Here a small DHPR trigger is more effective at high voltages (Figure 11D, dashed line).

The results demonstrate that this simplified model of functional release unit with Markov descriptions of DHPR and RyR channels produces both positive gain and gradedness when implemented as Monte Carlo simulations. The gradedness is generally robust with variation in the model parameters, and model behavior changes in reasonable and understandable ways (i.e., increases of SR load increase SR Ca²⁺ release). The model also gives some insight into the following areas: 1) the mechanism of release, 2) the role of SR Ca²⁺ depletion on termination of release, 3) the role of RyR adaptation on termination of release, and 4) possible physiological implications of geometry changes.

The results illustrate that RyR Ca²⁺ release is most effectively triggered by DHPR Ca²⁺ influx, but can also be triggered by Ca²⁺ release from neighboring RyRs (Fig. 10). The latter effect of self-regenerating RyR Ca²⁺ release is shown to decrease with lower SR loads or higher RyR adaptation rates (Figure 7 and Figure 8 and Figure 9 and Figure 10). The moderating effects of decreasing SR Ca²⁺ load and RyR adaptation work synergistically to keep SR unloading to ∼50%. RyR Ca²⁺ release is invariant with DHPR open time (Figure 10C).

The magnitude of the gain or amplification factor for CICR is similar to those measured experimentally (Wier et al). However, the shape of the gain curve produced by the model (Figure 5BD) differs from that observed experimentally (Wier et al). This deficiency of the model arises from the assumption that all RyR in the functional unit see the same [Ca²⁺]_SS. If spatial gradients of Ca²⁺ were present, higher-amplitude DHPR fluxes seen at lower depolarizations would recruit more RyRs through a greater spatial spread of Ca²⁺ to nearby RyRs than would the lower-amplitude DHPR fluxes seen at higher depolarizations. This would result in increased gain at lower depolarizations.

Previous modeling by our group that local depletion of SR Ca²⁺ may play an important role in terminating Ca²⁺ release. There are two proposed mechanisms. The first mechanism assumes the existence of two separate SR pools, the Ca²⁺ uptake pool (called network SR or NSR) and the Ca²⁺ release pool (called junctional SR or JSR). If Ca²⁺ transfer between these compartments is slow, then JSR can deplete, thus terminating RyR Ca²⁺ release (see (Jafri et al). However, the two-pool SR is a hypothetical construct, and there is little experimental support for such a long time constant between Ca²⁺ uptake and release sites (Bers, 1991). A second mechanism assumes that the NSR-JSR transfer rate is large, so that the [Ca²⁺] in each are very similar (i.e., SR is essentially one compartment). In this scheme, depletion occurs in both NSR and JSR (see (Jafri et al). A problem with this mechanism is that our previous simulations using a deterministic model suggest that almost complete depletion occurs during Ca²⁺ release. In contrast, experimental evidence suggests a maximal SR depletion of only ∼50% (Janczewski et al,Bassani et al).

The stochastic simulations here suggest that Ca²⁺ release termination can be achieved with a degree of total JSR unloading similar to that measured experimentally. In this model, the JSR is refilled by the NSR that has a fixed [Ca²⁺]. This is clearly a simplification; however, such a construction makes it clear that SR Ca²⁺ release can terminate in the absence of a large degree of total SR Ca²⁺ depletion in SR. Indeed, Ca²⁺ release is terminated in our model despite there being a continual refilling of JSR.

The mechanism of model Ca²⁺ release termination is the loss of self-regenerating Ca²⁺ release by the RyR as the local SR empties. In the numerical experiments where a single RyR is assumed to open, a large value of [Ca]_JSR produced a high degree of coupling between RyR openings and a correspondingly large Ca²⁺ release (Figure 10BD). In contrast, when [Ca]_JSR is reduced to 400μM, SR Ca²⁺ release is substantially diminished (Figure 10D). Hence the likelihood of self-regenerating RyR openings decreases as SR empties.

Without regeneration, the SR Ca²⁺ release can more closely follow DHPR openings, i.e., SR Ca²⁺ release occurs with DHPR openings and ceases shortly afterward. DHPR influx drives SR Ca²⁺ release for the complete duration of a voltage pulse, even after SR has depleted (i.e., see Figure 4B). From the single DHPR opening data, the SR Ca²⁺ release is approximately proportional to the SR load (Figure 10C). Therefore, at low SR loads, DHPR influx can effectively drive RyR Ca²⁺ release.

A proposed role for adaptation in the cardiac myocyte is to provide negative feedback to SR Ca²⁺ release to counter the strong positive feedback of CICR (Valdivia et al). We found that modifying adaptation rate did effect the magnitude of the Ca²⁺ release event, the temporal features of the Ca²⁺ transient, and the degree of unloading of SR (Figure 7B). However, at high SR Ca²⁺ load, release is sustained as a result of regenerative RyR opening despite the presence of adaptation (Figure 8B). Moreover, even if the adaptation rate is increased 10-fold, release does not terminate in the case of high SR Ca²⁺ load (data not shown). In fact, experiments show that under the conditions of high SR Ca²⁺ load, the mean open time of the RyR increases to as high as 17.4ms (Lukyanenko et al). They attribute this to modulation of RyR inactivation by SR [Ca²⁺]. In the simulations, doubling SR [Ca²⁺] increases the number of channel openings (from 14,520 to 29,490 for 500 runs) and the mean open time increases (from ∼2.8ms to ∼7.1ms) without the inclusion of any explicit effects of SR lumenal [Ca²⁺] on the RyR adaptation or inactivation. The mechanism is due to the differences in [Ca²⁺] at the mouth of the channel acting on the RyR and not on lumenal SR Ca²⁺ acting on RyR as observed by Xu and Meissner, 1998. Hence, we can conclude that SR emptying plays a more crucial role in terminating Ca²⁺ release than does adaptation. We note that others stochastic modeling of diadic SR release also suggest an important role of SR depletion in the termination of release (Stern, 1992).

Adaptation may also play other important roles in shaping RyR Ca²⁺ release. The data of Figure 10B show that integrated RyR flux is essentially independent of the duration of the DHPR open time. This lack of dependence is a consequence of RyR adaptation (the RyR moving from state P_O4 to state P_C5) so that longer duration triggers do not produce additional RyR Ca²⁺. This behavior is consistent with experimental findings that Ca²⁺ spark amplitudes are independent of the duration of the DHPR trigger current (Sham et al,Cannell et al) (see further discussion of Ca²⁺ sparks below). Other potential roles for adaptation could be to prevent secondary SR Ca²⁺ release events from occurring during an AP. Secondary SR Ca²⁺ release events are one potential mechanism for early afterdepolarizations, a proarrhythmogenic condition thought to be the basis of torsades de pointe and ectopic beating (Volders et al). Our previous work also suggests that the slow recovery from adaptation may play a crucial role in shaping mechanical restitution and other interval-force relations (Rice et al,Rice et al). More recently, it has been shown that mechanical restitution results from slow recovery of RyRs from the refractory state (Sham et al), a mechanism that has been suspected for quite some time (Bers, 1991).

Although not the primary focus of this paper, these simulations do bear on issues regarding properties of Ca²⁺ sparks. Because our simulations do not consider the spatial aspects of the functional unit, we cannot compare our results directly to experimental spark data. However, the magnitudes and temporal aspects of unitary Ca²⁺ release can be compared to experimental data.

The simulations that most closely match Ca²⁺ sparks are the RyR Ca²⁺ release events that are triggered by a single channel opening (Fig. 10). The size and shape of the [Ca²⁺] transients are similar to those reported in the literature (Sham et al). The half-time of decay for spark is measured to be 10 to 40ms (Santana et al). Similar to experimental data (Sham et al,Cannell et al), we find that Ca²⁺ spark amplitudes are independent of the duration of the DHPR trigger current, and the termination of SR Ca²⁺ releases does not depend on cessation of DHPR influx. The model is also consistent with the experimental observation that Ca²⁺ spark amplitudes are correlated with SR Ca²⁺ load (Satoh et al).

In our voltage pulse simulations, the initial Ca²⁺ release events from SR are large-amplitude, and thus are more likely the summation of multiple sparks. While the larger initial Ca²⁺ release appears to be regenerative, RyR Ca²⁺ release stops in response to the termination of DHPR current when the voltage step repolarizes (under default conditions). This is consistent with findings that halting DHPR currents can stop SR Ca²⁺ release (Wier et al,Cleemann and Morad, 1991).

In experiments, Gomez and co-workers observed that congestive heart failure decreases coupling between trigger DHPR Ca²⁺ influx and RyR Ca²⁺ release, possibly as a result of altered geometry in diadic space. (Gomez et al). Their results showed a bell-shaped curve with decreased gain across the voltage range tested. Similar results were not obtained in our simulations of altered diadic space geometry by either of two methods, increasing FRU subspace volume or increasing the efflux rate of Ca²⁺ out of the FRU subspace. In both cases, the RyR release was decreased at high voltages, consistent with the experimental finding, but the RyR release was increased at lower voltages, an effect not observed experimentally. Because neither of these two manipulations, either individually or together, reproduces the experimental data, the simulation results suggest that other changes besides these simple geometric alterations may be needed to account for the decreased gain observed in CHF. Indeed, the simulation results most similar to the experimental finds are for increased adaptation rate (Fig. 7, dashed trace) or for decreased SR load (Fig. 8, dotted trace). Eisner and co-workers, 1998 suggest that altering either the trigger flux from DHPR or the RyR flux alone will only have transient effects on SR Ca²⁺ release and that the only net result of reduced coupling would be a decrease in SR Ca²⁺ load. In fact, a decreased SR Ca²⁺ load is observed during heart failure, due to both down-regulation of the SERCA2a pump and up-regulation of Na⁺/Ca²⁺ exchange (O’Rourke et al,Winslow et al). However, Gomez and co-workers used a pulse protocol designed to control SR Ca²⁺ load in the experiments. CHF-induced changes in adaptation rate were not tested experimentally, and hence could be a potential mechanism of altered DHPR-RyR coupling during CHF.

While the model produces robust graded Ca²⁺ release, we cannot rule out other potential mechanisms that might also contribute to graded Ca²⁺ release in real cells. For example, Ca²⁺ diffusion in the diadic space might result in the recruitment of additional RyR in response to larger DHPR trigger Ca²⁺ influx. This mechanism could be termed recruitment within a functional unit as a result of local [Ca²⁺] gradients in microscopic domain. Another type of recruitment could also occur if cross-coupling exists between functional units. For example, a large Ca²⁺ release in one functional unit could potentially raise the [Ca²⁺] gradient high enough to produce Ca²⁺ release in a neighboring functional unit. There is some experimental evidence that Ca²⁺ release sites can be coherent over distances of 600nm that are larger than an FRU (Parker et al).

While we do assume a uniform [Ca²⁺] in each subspace, the need to consider each individual space is clear. The failure of the deterministic simulations here and elsewhere to produce graded Ca²⁺ release illustrates this point. This finding is predicted by previous work by Stern suggesting that “common pool” models, like our deterministic model, cannot produce both high gain and graded Ca²⁺ release without unrealistically tight control over parameters. In contrast, a “local control” model, like our stochastic model, can potentially reproduce graded Ca²⁺ release. Our results differ from previous work in that we show that robust high gain and graded Ca²⁺ release can occur in a model with detailed descriptions of DHPR and RyRs. Hence, these descriptions may be sufficient to explain graded Ca²⁺ release, although we cannot rule out other possible effects such as [Ca²⁺] gradients within a functional unit or cross-coupling between functional units.

A model is developed to describe the CICR in rat cardiac muscle. Using stochastic methods, CICR shows both high gain and gradedness similar to the experimental studies. The high gain and gradedness are generally robust when parameters are varied, although the amplitude and temporal details of Ca²⁺ release do change. The model suggest that DHPR influx produces a short, high-amplitude change in subspace [Ca²⁺] that is very effective in opening RyRs. For a single DHPR opening, the resulting RyR Ca²⁺ release flux is insensitive to the DHPR open duration but is strongly correlated with SR Ca²⁺ load, consistent with experimental Ca²⁺ spark data. In contrast, the single RyR openings require a long open duration and large SR load to be effective in triggering opening of neighboring RyRs. This low coupling between adjacent RyRs, especially as SR depletes, may explain why CICR does self-regenerate until the SR empties. Our results suggest that adaptation alone does not terminate Ca²⁺ release, but that adaptation plays an important modulatory role in shaping Ca²⁺ release duration and magnitude.

In contrast to the results above for stochastic simulations, our previous deterministic simulations show regenerative Ca²⁺ release only at optimum voltages for trigger influx (near 0mV), but no Ca²⁺ release outside this range (i.e., all-or-none response). The failure of our deterministic model suggests that CICR cannot be predicted by assuming a single subspace to represent the ensemble average of such spaces (“common pool model”). However, our stochastic results show that robust high gain and graded Ca²⁺ release are produced when the local [Ca²⁺] in each subspace is considered. Hence, our modeling results suggest that gradedness can occur by the recruitment of the independent functional units. Moreover, this gradedness occurs in the absence of spatial [Ca²⁺] gradients across the functional unit or cross-coupling between functional units, although our modeling cannot rule out these other potential mechanisms.