INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

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., 1998). 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., 1994; Lopez-Lopez et al., 1995). 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., 1993). We assume that single voltage-activated DHPR channels can open nearby RyRs via localized increases of [Ca²⁺] (Cannell et al., 1995; Lopez-Lopez et al., 1995). 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.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

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., 1998; Cannell et al., 1995). 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).

View larger version (20K): [in this window] [in a new window]   FIGURE 1   Schematic diagram of the functional release unit model. Ca2+ influx through the DHPR (JDHPR) enters the functional unit subspace (dashed box). This rise in local [Ca2+] fosters the opening of one or more of the eight RyR, producing additional influx (JRyR) 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 (Jxfer) occurs via simple diffusion, as does the flux to refill junctional SR (Jtr).

We also assume that each FRU behaves independently of others based on experimental findings (Sham et al., 1998; Cannell et al., 1995; Cheng et al., 1996). 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)

The Ca²⁺ flux (J_tr) that refills the junctional SR (JSR) is computed as the [Ca²⁺] gradient divided by the relaxation time constant (_tr):

(2)

where [Ca²⁺]_JSR is the JSR Ca²⁺ concentration and [Ca²⁺]_NSR is the network SR (NSR) Ca²⁺ concentration. It is assumed that the activation of an individual FRU does not have a significant impact on the total [Ca²⁺]_NSR. Hence, [Ca²⁺]_NSR is a constant ([Ca²⁺]_NSR = 800 µM in the default model).

The DHPR channel is represented by a mode-switching model that was developed previously (Jafri et al., 1998). A state diagram of the mode switching and voltage activation is shown in Fig. 2 A. 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.

View larger version (26K): [in this window] [in a new window]   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 Cn, where n is the number of permissive subunits. With depolarization, the channels undergo transitions from left to right. With four permissive subunits (C4), there is a voltage-independent transition to the conducting state O. With elevation of [Ca2+] the channel transition occurs from the Mode Normal (top) to Mode Ca (bottom). f f' so that transitions into OCa are rare, which effectively inactivates the channel in Mode Ca. (B) A state model for the RyR with adaptation. The states PC1 and PC2 are the resting states, PO3 and PO4 are the open states, and PC5 and PC6 are the adapted or refractory states. The transitions governed by the rates k1,2, k2,3, k3,4, k5,4, k5,6, and k2,5 are all Ca2+-dependent (see Table 2 for values of the rate constants).

Transitions to Mode Ca are controlled by , which is a function of Ca²⁺. Moving rightward in Fig. 2 A, 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 I. 

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., 1999). 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 Fig. 2 B. 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.

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

(5)

where _RyR 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 (Wagner and Keizer, 1994) with _SS describing buffering in the subspace and _JSR describing buffering in the JSR. In the subspace there is buffering by the negatively charged phospholipid headgroups in the sarcolemma and the SR membrane. These are approximated as having fixed concentrations, and Ca²⁺ dissociation constants are given by [B]_SL and K_BSL for the sarcolemma and [B]_SR and K_BSR for the SR (Smith et al., 1998). In the JSR, calsequestrin (CSQN) is assumed to be the only significant Ca²⁺ buffer (Luo and Rudy, 1994).

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.1 s) to reach steady state at a holding potential of 80 mV, followed by a voltage step to a test potential for 0.2 s, and then returned to the 80 mV holding potential. Each trial is repeated 500 times before changing the voltage step (typically increased from 50 and +50 mV 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.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

Sample Monte Carlo results are shown in Fig. 3 for a single trial with a voltage step to 0 mV. The peak fluxes are shown for the DHPR channel (solid line) and the sum of the eight RyRs (dashed line) in Fig. 3 A. After a voltage transition to 0 mV at time 0.1 s, the DHPR might open (Fig. 3 A). 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 Fig. 3 B, 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. Fig. 3, C and D 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.

View larger version (18K): [in this window] [in a new window]   FIGURE 3   Sample results for a single functional unit driven by a voltage step to 0 mV from 0.1 to 0.3 s. (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 Ca2+ concentration in the subspace ([Ca2+]SS). (D) The changes in JSR Ca2+ concentration ([Ca2+]JSR) indicates the degree to which SR empties during a Ca2+ 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 Fig. 4, A and B, 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. Fig. 4, C and D 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., 1995a).

View larger version (21K): [in this window] [in a new window]   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 0 mV from 0.1 to 0.3 s. (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 Ca2+ concentration in the subspace ([Ca2+]SS). (D) The changes in average JSR Ca2+ concentration ([Ca2+]JSR) indicate the degree to which SR empties during a Ca2+ 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 +50 mV in 5-mV increments. When plotted against clamp voltage, the averaged DHPR fluxes produce bell-shaped profiles for both integrated (Fig. 5 A, solid line) and peak fluxes (Fig. 5 C, solid line). The peak of the DHPR flux is ~10 mV, 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 Fig. 5 B and the gain for the peak flux is shown in Fig. 5 D. While the gain shown is in the range observed physiologically (Wier et al., 1994), the shape differs from that seen in the experiments. Possible explanations for this will be described in the Discussion.

View larger version (24K): [in this window] [in a new window]   FIGURE 5   Model exhibits positive gain and graded Ca2+ release. To modulate the amount of trigger Ca2+, voltage pulses are changed in 5-mV steps from 50 to 50 mV. 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 Ca2+ 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 Fig. 6 A, the activation function of the DHPR channel is shifted by 10 and +10 mV, producing shifts of the rising edge of the peak DHPR flux functions (seen by the three low-amplitude curves in Fig. 6 A). The rising phases of the RyR peak flux functions show corresponding shifts (seen by the three higher amplitude curves in Fig. 6 A). 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 Fig. 6 B, shifting the open-channel I-V relation by 10 and +10 mV 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.

View larger version (40K): [in this window] [in a new window]   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 Ca2+, voltage pulses are changed in 5-mV steps from 50 to 50 mV. The data are averaged from 500 runs. (A) The activation function of the DHPR channel is shifted by 10 (dashed line) and +10 mV (dotted line). The rising phase of the RyR average peak flux curves shows corresponding shifts, showing that RyR Ca2+ 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 +10 mV (dotted line). The falling phase of the RyR average peak flux curves show corresponding shifts, showing that RyR Ca2+ 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 Fig. 2 B). A slower adaptation rate increases peak RyR flux and [Ca²⁺] in the subspace (Fig. 7, A and B, dotted lines). Conversely, the effect of more rapid adaptation is to decrease both peak flux and [Ca²⁺]_SS (Fig. 7, A and B, 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., 1994). This finding suggests that measurements of gain may help to limit estimates of RyR adaptation rates.

View larger version (27K): [in this window] [in a new window]   FIGURE 7   The effects of modifying adaptation rate on RyR Ca2+ release. The transition rates into the adapted state (k2,5 and k4,5 in Fig. 2 B) 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 50 mV 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 Fig. 5 C). Increased adaptation rate decreases the peak flux. (B) The average [Ca2+]SS transients at 0 mV show that adaptation rate effects both the magnitude and temporal profile of RyR Ca2+ release. (C) The average [Ca2+]JSR lines show that the higher rate of adaptation (dashed line) decreases the degree of SR unloading during Ca2+ release, whereas decreased rate (dotted line) produces greater SR unloading.

The effect of RyR adaptation on JSR Ca²⁺ concentration is examined in Fig. 7 C. 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., 1995a; Janczewski et al., 1995). For example, by decreasing the SR Ca²⁺ content by 56%, its control value reduces the RyR Ca²⁺ release by 52% (Janczewski et al., 1995). 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 (Fig. 8 A, dashed line). Similarly, when SR Ca²⁺ load is halved from the control value (400 µM vs. 800 µM), the peak flux decreases (Fig. 8 A, dashed line).

View larger version (27K): [in this window] [in a new window]   FIGURE 8   The effects of modifying SR Ca2+ load on RyR Ca2+ release. The [Ca2+] of JSR and the [Ca2+] 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 50 mV 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 Fig. 5 C). (B) The average [Ca2+]SS transients at 0 mV show that SR load effects both the magnitude and temporal profile of RyR Ca2+ release. Doubling the SR load (dashed line) produces elevated sustained Ca2+ release during the voltage pulse (0.1-0.3 s) and regenerative RyR Ca2+ release after the pulse when DHPR influx has ceased. (C) The average [Ca2+]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 Fig. 8 B, 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.3 s in Fig. 8 B). 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 Fig. 8 B 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 Figs. 3-5. In Fig. 9 A, the DHPR opens times are plotted during the 0.1-0.3-ms voltage step to 0 mV. Most points are below 3 ms, 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 Fig. 9 A. The line has slope near zero, demonstrating that the DHPR mean open time is constant at ~0.5 ms. This corresponds well with the value of 0.45 ms measured experimentally by Rose and co-workers (1992). In Fig. 9 B, the RyR open times show some time dependence with a peak mean open time of ~6.7 ms early during the voltage pulse and mean open time of ~2.8 ms later in the pulse. These values are in agreement with the ranges of values from 1.5 to 7.1 ms measured experimentally for control conditions (Lukyanenko et al., 1996; Eager and Dulhunty, 1998). For 500 runs, the 14,520 RyR openings (Fig. 9 B) are much more numerous than 1567 DHPR openings (Fig. 9 A).

View larger version (29K): [in this window] [in a new window]   FIGURE 9   The DHPR and RyR open times during a voltage pulse to 0 mV 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.51 ms for 1567 points collected from 500 runs. (B) RyR open times show some time-dependence with a peak average open time of ~6.7 ms early during the voltage pulse, and average around 2.8 ms (white line) later in the clamp. The mean open time is 2.8 ± 3.8 ms 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 Fig. 10 A for a DHPR open time of 0.5 ms (black line) and 16.0 ms (gray line). The amplitudes are similar with a slight variation of the time course. This compares well with the experimental results of Sham and co-workers (1998) in which the mean open time of the DHPR is increased from 0.27 ms to 15.9 ms 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.5 ms (black line) and 16.0 ms (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 (Figs. 3-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.6 s, data not shown). This finding suggests that SR unloading plays an important role in the termination of SR Ca²⁺ release.

View larger version (28K): [in this window] [in a new window]   FIGURE 10   The RyR Ca2+ release elicited from a single DHPR or RyR opening. (A) The average Ca2+ transients for 200 runs with default initial SR Ca2+ load ([Ca2+]JSR = 800 µM). The response for a 0.5 ms (black line) and a 16 ms (gray line) DHPR opening are similar, although the latter is slightly larger and has a small shoulder as a result of the DHPR Ca2+ influx. (B) The Ca2+ transients for a single a 0.5 ms (black line) or 16 ms (gray line) RyR opening are similar, although slightly smaller that those produced by a single DHPR opening. (C) RyR Ca2+ release is independent of the open duration of a single DHPR, but highly dependent on SR load. The [Ca2+] of JSR and the [Ca2+] of the refilling store are set to 800, 600, 400, and 200 µM (as labeled) to vary the SR load. (D) RyR Ca2+ release shows some dependence on the open duration of a single RyR for short open times. As in (C) above, RyR Ca2+ 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 Fig. 4 A, 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 (Fig. 9 B), 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 (Fig. 10 B). 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 Fig. 10, C and D. The peak RyR Ca²⁺ flux for different values of initial [Ca]_JSR are shown. In Fig. 10 C, the DHPR was held open for a range of open times from 0.1 ms to 50.0 ms. The peak RyR Ca²⁺ flux is invariant of DHPR open time, as observed by Sham and co-workers (1998), 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 Fig. 10 D, 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., 1997). 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 (Fig. 11 A, dashed line) produces a larger peak flux than control conditions (solid line, same as Fig. 5 C), 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 (Fig. 11 A, 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 (Fig. 11 B, dotted line), and the depletion of Ca²⁺ from the SR is less for the smaller subspace volume (not shown).

View larger version (30K): [in this window] [in a new window]   FIGURE 11   The effects of modifying subspace volume on RyR Ca2+ release. The volume of the subspace is doubled (dashed line) or halved (dotted line). Data are from voltage pulses from 50 to 50 mV 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 Fig. 5 C), but only at lower voltages. At higher voltages, the trigger influx through the DHPR channel is insufficient to trigger Ca2+ release. The smaller subspace volume (dotted line) produces greater Ca2+ 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 [Ca2+]SS transients at 0 mV show that subspace volume effects both the magnitude and temporal profile of RyR Ca2+ release. Halving the subspace (dotted line) produces elevated sustained Ca2+ release during the voltage pulse (0.1 to 0.3 s) and a small amount regenerative RyR Ca2+ release after the pulse when DHPR influx has ceased. Note that doubling the subspace increased peak flux at 0 mV (A), but produced a smaller change in [Ca2+]SS given the larger volume. (C) Doubling the transfer rate out of the subspace (xfer) results decreases RyR Ca2+ release at high depolarizations (dotted line). Decreasing xfer by half increases the gain of RyR Ca2+ at higher voltages. (D) The average [Ca2+]SS transients at 0 mV 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 (Fig. 11 C, dotted line). Similar to the increased subspace volume case, the doubled transfer rate decreases the effectiveness of the DHPR trigger at high voltages (Fig. 11 D, dotted line). Reducing the transfer rate by half causes the right side of the graded release curve to be shifted to the right (Fig. 11 C, dashed line). Here a small DHPR trigger is more effective at high voltages (Fig. 11 D, dashed line).

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

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.

Mechanism of release

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 (Figs. 7, 8, and 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 (Fig. 10 C).

The magnitude of the gain or amplification factor for CICR is similar to those measured experimentally (Wier et al., 1994). However, the shape of the gain curve produced by the model (Fig. 5, B and D) differs from that observed experimentally (Wier et al., 1994). 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.

Depletion of SR and termination of Ca²⁺ release

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., 1998). 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., 1998). 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., 1995; Bassani et al., 1995b).

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 (Fig. 10, B and D). In contrast, when [Ca]_JSR is reduced to 400 µM, SR Ca²⁺ release is substantially diminished (Fig. 10 D). 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 Fig. 4 B). From the single DHPR opening data, the SR Ca²⁺ release is approximately proportional to the SR load (Fig. 10 C). Therefore, at low SR loads, DHPR influx can effectively drive RyR Ca²⁺ release.

Role of adaptation and the termination of 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., 1995). 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 (Fig. 7 B). However, at high SR Ca²⁺ load, release is sustained as a result of regenerative RyR opening despite the presence of adaptation (Fig. 8 B). 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.4 ms (Lukyanenko et al., 1996). 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.8 ms to ~7.1 ms) 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 Fig. 10 B 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., 1998; Cannell et al., 1995) (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., 1997). 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., 1998, 1999). More recently, it has been shown that mechanical restitution results from slow recovery of RyRs from the refractory state (Sham et al., 1999), a mechanism that has been suspected for quite some time (Bers, 1991).

Relationship to Ca²⁺ sparks

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., 1998). The half-time of decay for spark is measured to be 10 to 40 ms (Santana et al., 1996). Similar to experimental data (Sham et al., 1998; Cannell et al., 1995), 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., 1997).

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., 1994; Cleemann and Morad, 1991).

Physiological implications of geometry changes

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., 1997). 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., 1999; Winslow et al., 1999). 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.

Other potential mechanisms of graded Ca²⁺ release

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 600 nm that are larger than an FRU (Parker et al., 1996).

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.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUSIONS REFERENCES

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 0 mV), 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.