INTRODUCTION

Top Abstract Introduction Methods Results Concluding Remarks Appendix References

There are several models for cardiac action potentials (APs) that include intracellular Ca²⁺ handling (DiFrancesco and Noble, 1985; Luo and Rudy, 1994a, b; Nordin, 1993; Lindblad et al., 1996). The two most well known are the Luo-Rudy Phase II ventricular cell model and the DiFrancesco-Noble Purkinje cell model. Whereas these models generate APs using biophysically detailed descriptions of membrane currents, the calcium subsystem is represented by a phenomenological model that mimics Ca²⁺-induced Ca²⁺ release (CICR), but fails to capture the biophysical detail of the mechanisms involved. Other models exist that describe only the Ca²⁺ handling aspects of the cardiac myocyte (Dupont et al., 1996; Tang and Othmer, 1994; Schouten et al., 1987). However, these models do not include the membrane currents, and thus, cannot capture the interplay between the Ca²⁺ subsystems and membrane currents. Such shortcomings have limited the predictive ability of the models with regard to Ca²⁺ dynamics and force generation.

New experimental evidence concerning the RyR and L-type Ca²⁺ current now provides a basis for biophysical models of Ca²⁺ handling in cardiac cells, and for the integration of these models into existing descriptions of the cardiac AP. Recently, Györke and Fill (1993) used lipid bilayer preparations to show that the RyR displays "adaptation" in response to successive Ca²⁺ elevations. Adaptation means that incremental release of Ca²⁺ from RyR occurs in response to incremental step changes in cytoplasmic side Ca²⁺. Mechanisms of adaptation have been postulated and presented in mathematical models by several groups (Keizer and Levine, 1996; Sachs et al., 1995; Cheng et al., 1995; Tang and Othmer, 1994).

Clearly, properties of CICR and adaptation of RyR will have an important influence on the detailed shape of the cardiac cell Ca²⁺ transient. The shape of this Ca²⁺ transient is also influenced by properties of Ca²⁺-mediated inactivation of the L-type Ca²⁺ channel (Grantham and Cannell, 1996). Recently, Imredy and Yue (1994) have provided direct evidence concerning the molecular mechanism of this inactivation. Specifically, they have proposed that as Ca²⁺ binds to the L-type Ca²⁺ channel, the channel undergoes a mode switch from one characterized by dense bouts of activity (mode normal) to one typified by infrequent openings (mode Ca). This results in time-dependent inactivation kinetics that are very different from those seen in previous cardiac Ca²⁺ channel models (DiFrancesco and Noble, 1985; Luo and Rudy, 1994a, b; Nordin, 1993).

In this paper we formulate a multimodal model of Ca²⁺-mediated inactivation of the L-type Ca²⁺ channel based on the concept of mode switching. We incorporate this model, along with a modified version of the Keizer-Levine RyR adaptation model and membrane currents from the Luo-Rudy Phase II model, into a model for ventricular cell Ca²⁺ handling. We then use the model to explore the effects of RyR adaptation and SR load on the cardiac action potential. We also investigate the interplay of RyR adaptation and SR load on the frequency-dependent aspects of the cardiac Ca²⁺ transient.

	MECHANISMS

Our model is based on the Luo-Rudy Phase II model for ventricular action potentials with several modifications (Fig. 1): 1) the Luo-Rudy Phase II L-type Ca²⁺ current is replaced with our new formulation based on the mode-switching behavior observed in rats by Imredy and Yue (1994); 2) the Luo-Rudy Phase II Ca²⁺ SR release mechanism is replaced with the Keizer-Levine RyR model with adaptation (Keizer and Levine, 1996) based on data from isolated canine RyRs; 3) the RyRs and L-type Ca²⁺ channels are assumed to empty into a restricted subspace located between the junctional sarcoplasmic reticulum (JSR) and T-tubules; 4) both high- and low-affinity Ca²⁺ binding sites for troponin are included; 5) the magnitude of the Luo-Rudy Phase II membrane currents I_Kp, I_Na, I_NaCa, I_ns(Ca), and I_Ca,b are scaled to preserve myoplasmic ionic concentrations and AP shape. Table 6 (Appendix 3) shows the parameters for the membrane currents. Those that were rescaled have their original Luo-Rudy Phase II values shown in parentheses.

View larger version (44K): [in this window] [in a new window]   FIGURE 1   Schematic diagram of mechanisms involved in cardiac Ca2+ dynamics. In response to membrane depolarization during an AP, L-type Ca2+ channels open, allowing the influx of Ca2+ (ICa) into a restricted subspace, where it triggers Ca2+ release from the JSR (Jrel) via ryanodine receptors (RyR). Ca2+ diffuses from the subspace to the bulk myoplasm (Jxfer), where it is removed from the cell by Na+-Ca2+ exchangers (INaCa) and ATP-dependent Ca2+ pumps (Ip(Ca)) or resequestered into NSR by SERCA pumps (Jup). The JSR is refilled by diffusion from the NSR (Jtr). Ca2+ is buffered by calmodulin in the subspace and myoplasm, by troponin in the myoplasm, and by calsequestrin in the JSR. There is a leak (Jleak) from the NSR and a background Ca2+ current across the sarcolemma (included in Imembrane), which help to maintain Ca2+ homeostasis. Imembrane includes all of the other transarcolemmal currents that are not defined above.

L-type Ca²⁺ channel

Ca²⁺-mediated inactivation of the L-type channels has been shown recently to result from an intrinsic channel Ca²⁺-binding motif called an EF-hand (de Leon et al., 1995). The discovery of this Ca²⁺-binding site is important because it suggests that inactivation depends on local Ca²⁺ concentration instead of other hypothesized extrinsic mechanisms such as phosphorylation or second messengers. Inactivation occurs as Ca²⁺ binding induces the channel to switch (from mode normal) to a mode in which transitions to open states are extremely slow (mode Ca), as described in a model by Imredy and Yue (1994). It should be noted that this method of inactivation differs from that of current cardiac cell models. In these models, inactivation is a saturating Michaelis-Menten-type function of [Ca²⁺]_i that can instantaneously prevent channel permeation. In contrast, the model of Imredy and Yue predicts that inactivation via mode switching occurs with some delay after channel activation. Similarly, the model predicts that recovery to mode normal is not instantaneous when Ca²⁺ concentration falls.

The data of Imredy and Yue show that with an increase of local Ca²⁺, the L-type channel shifts to a gating mode that shows very infrequent openings. The original Imredy-Yue model described this two-mode behavior with five states (three in mode normal and two in mode Ca), and serves as the starting point of an improved L-type channel model described here. To complete the original Imredy-Yue channel model, the following features are added: 1) additional states; 2) improved Ca²⁺ inactivation, 3) voltage-dependent activation; 4) voltage-dependent inactivation; and 5) open-channel ion permeation.

A state diagram for mode switching and voltage-dependent activation is shown in Fig. 2. The upper and lower rows of states comprise the 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 or reverse rates between the closed states is dictated by the four-way symmetry assumed for the channel subunits. Voltage-dependent activation is incorporated through  and , which are increasing and decreasing functions of voltage. 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 (14K): [in this window] [in a new window]   FIGURE 2   Schematic diagram for the transitions between the states of the L-type Ca2+ channel. The upper row of states comprises the mode normal, and the lower row comprises mode Ca. The channel is composed of 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, there is a voltage-independent transition to the conducting state O. With elevation of [Ca2+]ss 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.

The transitions to mode Ca are controlled by , which is a function of Ca²⁺. As one moves right in Fig. 2, 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 (Shirokov et al., 1993; Hadley and Lederer, 1991). Additional constraints on the rate constants come from thermodynamic microreversibility. Microreversibility requires that for each cycle, the product of rates is equal whether taken in the clockwise or the counterclockwise direction. A similar formulation and state diagram has been proposed for N-type Ca²⁺ channels in bullfrog sympathetic neuron (Cory et al., 1993).

Voltage-dependent inactivation is modeled as a Hodgkin-Huxley-type gate. This gate can inactivate the channel independently of the states already discussed above. The time constant of inactivation is significantly longer than that of existing models of L-type current (DiFrancesco and Noble, 1985; Luo and Rudy, 1994a). These previous models probably underestimated the time constant of voltage-dependent inactivation, because the temporal properties of Ca²⁺-induced inactivation were not appropriately considered. That is, Ca²⁺-induced inactivation was assumed to be an instantaneous function of [Ca²⁺]_i with dynamics due solely to voltage-gated inactivation.

L-type channels have complex permeation properties, showing extremely high selectivity coupled with high open channel flux rates (Hess, 1990; Hille, 1992). To account for these properties, researchers have proposed a number of multiple binding site models (Hess and Tsien, 1984). However, these permeation models were not chosen for the improved Ca²⁺ subsystem because 1) their relatively high complexity would substantially increase computation and data storage demand; and 2) none of the models are consistent with all of the experimental data (Hess, 1990). A much simpler approach, used in existing cardiac models (Luo and Rudy, 1994a; DiFrancesco and Noble, 1985; Nordin and Ming, 1995), makes use of a Goldman-Hodgkin-Katz (GHK) formalism with large permeability for Ca²⁺ coupled to smaller permeabilities for K⁺ and Na⁺. However, this approach is also unsatisfactory because 1) the open channel current-voltage (I-V) relations do not match experimental results; and 2) GHK assumes that ionic species permeate independently (Hille, 1992), whereas data indicate essentially no monovalent permeation while the channel is passing significant amounts of Ca²⁺ current. As a compromise between the complex multisite models and the simpler GHK formalism, a novel permeation model is employed. In this model, it is assumed that 1) the Ca²⁺ current follows constant field theory and is the only inward current passing through the channel; and b) the permeability of K⁺ is a decreasing function of Ca²⁺ current. Under these assumptions, it becomes increasingly hard for K⁺ to pass when the channel is occupied by Ca²⁺ ions. The P_K' (the permeability of K⁺ as modified by the Ca²⁺ current I_Ca) is given by

(1)

where _K is the permeability of K⁺ in the absence of Ca²⁺ current, I_Ca is L-type Ca²⁺ current, and I_Cahalf is the level of Ca²⁺ current that reduces the permeability of K⁺ by 50%.

Open-channel I-V relations for this permeation model are shown by the solid trace in Fig. 3 A. The dashed traces show the separate contributions from the Ca²⁺ and K⁺ currents. In Fig. 3 B, the open-channel I-V relation is multiplied by the open probability (ignoring voltage and Ca²⁺-induced inactivation) to effectively simulate voltage clamp data, where peak current is measured before the channel inactivates significantly. The model data are similar to experimentally determined I-V relations (McDonald et al., 1986).

View larger version (21K): [in this window] [in a new window]   FIGURE 3   (A) Open L-type Ca2+ channel I-V relation for the permeation model (solid line), and the separate contributions by Ca2+ influx and K+ efflux (dashed lines). (B) Peak ICa obtained by multiplying the open channel I-V relation from above with the maximum open channel probability.

Ryanodine receptors

Ca²⁺ release from the JSR is based on a model for RyR by Keizer and Levine (Fig. 4). This model has two open states (P_O1 and P_O2) and two closed states (P_C1 and P_C2). At rest, the channel resides primarily in the first closed state, P_C1. Upon an increase in Ca²⁺, the channel switches briefly to the first open state P_O1, allowing Ca²⁺ to move through the channel, before it adapts by its transition to P_C2. Upon additional increases in Ca²⁺, the channel reopens by its transition to state P_O2, displaying the adaptive behavior seen experimentally (Györke and Fill, 1993). More recent experimental evidence (Valdivia et al., 1995) suggests that in the presence of physiological levels of cytosolic Mg²⁺, the rate of adaptation is faster than the findings of the original experiments by Györke and Fill, which lacked Mg²⁺. We follow the suggestion of Keizer and Levine that a faster rate for k_c⁺ would account for this finding. The bilayer experiments by Györke and Fill were performed at room temperature (25°C), whereas our model describes a cardiac ventricular cell at 37°C. To account for this difference, all of the rate constants are increased over the original Keizer and Levine values.

View larger version (13K): [in this window] [in a new window]   FIGURE 4   Schematic diagram of the transitions between the states of the RyR. PC1 and PC2 are closed states, and P01 and P02 are open states. The k's are the transition rates. The transitions from PC1 to P01 and from P01 to P02 are Ca2+ dependent. At resting Ca2+ levels, most of the channels are in state PC1. With an increase in Ca2+, the channels undergo transitions to P01 before channel adaptation by transition to PC2. Additional increases in Ca2+ cause reopening by transitions into P02.

The RyR is also modified to account for the environment of the subspace. The Keizer-Levine model assumes that the RyR can be exposed to peak [Ca²⁺]_i (Appendix 2) values of ~1.0 µM. In our model, the RyR is located in the subspace, where it is exposed to [Ca²⁺]_ss in excess of 10.0 µM (Santana et al., 1996; Cannell et al., 1994; Isenberg, 1995). In addition, [Ca²⁺]_ss levels change more rapidly than [Ca²⁺]_i. Thus the rate constants are modified to adjust the channel sensitivity to Ca²⁺, so that the channel functions properly in the appropriate Ca²⁺ range (Table 3, Appendix 3.3).

Subspace

Many existing cardiac AP models have Ca²⁺ influx from an L-type channel and release from SR via RyR emptying into the bulk myoplasm (DiFrancesco and Noble, 1985; Luo and Rudy, 1994a, b). In contrast, many researchers have suggested that both of these channels empty into a more restricted subspace, where Ca²⁺ concentration may increase to much greater levels than in the bulk myoplasm (Isenberg, 1995; Nordin, 1993; Stern, 1992; Bers, 1991; Lederer et al., 1990). The anatomical evidence suggests that L-type and RyR channels exist in a ratio of 1:5.6 in the guinea pig (Bers and Stiffel, 1992). The channels are located in specialized junctional areas, where T-tubular membrane and SR membrane are separated by a junctional gap of 12 nm (Frank, 1990). Following Isenberg (1995), we define a functional unit as four L-type channels surrounded by 25 RyRs. Assuming that the RyRs are arranged in 5 × 5 array according to a spacing of 60 nm, the cluster will cover a surface area of 300 nm × 300 nm. The number of functional units can be estimated as 5500 L-type channels per cell (Isenberg, 1995) divided by four channels per functional unit. The volume of the restricted subspace can be estimated by

(2)

The volume of the subspace is then calculated as 1.485 × 10⁹ µl. For comparison, this volume is about four orders of magnitude smaller than the total volume of the myoplasm (25.84 × 10⁶ µl).

The rate of Ca²⁺ flux out of this subspace (J_xfer) allows [Ca²⁺]_ss and [Ca²⁺]_i to equilibrate with a time constant _xfer (Eq. 8).

Buffering

Ca²⁺ buffers sequester in excess of 98% of the total Ca²⁺ released during a contraction (Berlin et al., 1994). These buffers are distributed between mobile and stationary buffers. In cardiac cells, the mobile buffer is calmodulin, and the stationary buffers are troponin and myosin (Robertson et al., 1981). In previous models, the effects of myosin have been ignored because, in simulations by Robertson and co-workers, it was never more than 2.3% saturated during a Ca²⁺ transient (Robertson et al., 1981). Calsequestrin is the most abundant SR luminal Ca²⁺-binding protein (Isenberg, 1995). Immunoelectron microscopy has shown that calsequestrin is located in the JSR and not in the NSR (Jorgensen and Campbell, 1984).

Buffers B bind Ca²⁺ by the following equation:

(3)

If the kinetics of binding is fast compared to Ca²⁺ release, Eq. 3 remains close to equilibrium during the Ca²⁺ transient. In this case, the rapid buffering approximation is used (Wagner and Keizer, 1994). Calsequestrin and calmodulin are considered to be fast buffers, based on their rate constants (Robertson et al., 1981; Cannell and Allen, 1984). Otherwise, Ca²⁺ binding to buffers lags behind the equilibrium values, as is the case with troponin. In this case, differential equations based on reaction kinetics must be solved. Both the high- and low-affinity sites for troponin are used in this model.

	MODEL EQUATIONS

Ca²⁺ subsystem

The Ca²⁺ subsystem consists of 10 differential equations specifying the time rate of change in 1) the myoplasmic Ca²⁺ concentration ([Ca²⁺]_i); 2) the subspace Ca²⁺ concentration ([Ca²⁺]_ss); 3) the junctional SR Ca²⁺ concentration ([Ca²⁺]_JSR); 4) the network SR Ca²⁺ concentration ([Ca²⁺]_NSR); 5) the probability that RyR is in the first closed state (P_C1); 6) the probability that RyR is in the first open state (P_O1); 7) the probability that RyR is in the second open state (P_O2); 8) the probability that RyR is in the second closed state (P_C2); 9) the concentration of Ca²⁺ bound to high-affinity troponin-binding sites ([HTRPNCa]); and 10) the concentration of Ca²⁺ bound to low-affinity troponin-binding sites ([LTRPNCa]). The other buffers, calsequestrin in the JSR, and calmodulin in the myoplasm and subspace ([CSQN]_JSR, [CMDN]_i, and [CMDN]_ss, respectively), are fast buffers, so the rapid buffering approximation is used to minimize computation (Wagner and Keizer, 1994).

There are six Ca²⁺ fluxes to consider in the subsystem: 1) the Ca²⁺-induced Ca²⁺ release (CICR) flux via the RyR (J_rel; Eq. 4); 2) the leak flux from the NSR to the myoplasm (J_leak; Eq. 5); 3) the Ca²⁺ uptake flux by the SR Ca²⁺-ATPase pump (J_up; Eq. 6); 4) the transfer flux of Ca²⁺ from the NSR to the JSR (J_tr; Eq. 7); 5) the transfer flux from the subspace to the myoplasm (J_xfer; Eq. 8); and 6) the buffering flux of Ca²⁺ binding to troponin (J_trpn; Eq. 9). Three additional membrane current fluxes are necessary for the formulation of Ca²⁺ regulation: 7) the pump current (I_p(Ca); Eq. 84) flux of Ca²⁺ removal from the cell; 8) the L-type Ca²⁺ channel Ca²⁺ flux (I_Ca; Eq. 37); and 9) the Na⁺-Ca²⁺ exchange current (I_NaCa; Eq. 75) Ca²⁺ flux.

The RyR channel flux is described by

(4)

where v₁ is the maximum RyR channel Ca²⁺ flux, and RyR_open is the sum of the fraction of channels in the RyR channel open states P_O1 and P_O2 (RyR_open = P_O1 + P_O2). The equations describing RyR are given in Appendix 2. The leak from the NSR into the myoplasm is given by

(5)

where v₂ is the Ca²⁺ leak rate constant from the NSR. Ca²⁺ uptake into the NSR by the SR Ca²⁺-ATPase pump is given by

(6)

where v₃ is the SR Ca²⁺-ATPase (SERCA2a) maximum pump rate, and K_m,up is the half-saturation constant for the SR Ca²⁺-ATPase pumps. In this model, the SR Ca²⁺-ATPase pump has a Hill coefficient of 2, consistent with recent experimental findings (Lytton et al., 1992). The transfer flux of Ca²⁺ from the JSR to the NSR is described by

(7)

where _tr is the time constant for transfer from NSR to JSR. The transfer flux from the subspace to the myoplasm is

(8)

where _xfer is the time constant for transfer from subspace to myoplasm. Buffering by troponin is described by

(9)

where [LTRPN]_tot is the total myoplasmic troponin low-affinity site concentration, [HTRPN]_tot is the total myoplasmic troponin high-affinity site concentration, k_htrpn⁺ is the Ca²⁺ on rate constant for troponin high-affinity sites, k_htrpn is the Ca²⁺ off rate constant for troponin high-affinity sites, k_ltrpn⁺ is the Ca²⁺ on rate constant for troponin low-affinity sites, and k_ltrpn is the Ca²⁺ on rate constant for troponin low-affinity sites. The first term in Eq. 9 describes the rate of Ca²⁺ binding to troponin high-affinity Ca²⁺-binding sites. The second term describes the off rate from these sites. The third and fourth terms (second line) describe the on and off rates for Ca²⁺ binding to troponin low-affinity binding sites, respectively.

The buffering effects of calmodulin and calsequestrin can be approximated by the rapid buffering approximation (Wagner and Keizer, 1994). This simplification assumes that the free and bound Ca²⁺ buffers with fast kinetics are at equilibrium. Using this assumption, a Ca²⁺-dependent function can be found that will scale the fluxes that comprise the time rate of change in total free Ca²⁺ by the factors B_i, B_ss, and B_JSR for the myoplasm, subspace, and JSR, respectively:

(10)

(11)

(12)

The balance equation for [Ca²⁺]_i is

(13)

where A_cap is the capacitive membrane area, V_myo is the myoplasmic volume, 2 is the valence of Ca²⁺, and F is Faraday's constant. Changes to [Ca²⁺]_ss are determined by

(14)

Note that the fluxes J_rel and J_xfer are scaled for the volume of the compartment. A detailed description of L-type Ca²⁺ current is presented in the next section. The [Ca²⁺]_JSR and [Ca²⁺]_NSR balance equations are

(15)

(16)

L-type Ca²⁺ channel

We created a new mathematical model to describe the L-type channel that is based on the experimentally observed mode-switching behavior of the channel. The rate constants for the L-type Ca²⁺ channel depend on the parameters in Table 4 (Appendix 3) and

(17)

(18)

(19)

(20)

(21)

where , , , ', and ' are in units of ms¹. This yields the following set of equations for channel states:

(22)

(23)

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

The state O in Eq. 27 is the conducting state. Equations 22-27 describe the normal mode, and Eqs. 28-33 describe mode Ca²⁺. The voltage inactivation gate y of the L-type Ca²⁺ channel is determined by

(34)

where

(35)

and

(36)

The L-type Ca²⁺ channel can then be written by

(37)

where _Ca is the maximum L-type Ca²⁺ channel conductance, y is the voltage inactivation gate, {O + O_Ca} is the open probability based on the mode-switching model (Eqs. 27 and 33), and 4 corresponds to the square of the charge of Ca²⁺. The last expression in Eq. 37 is the electrochemical driving force for L-type Ca²⁺ current. The first term in the numerator of this expression uses 0.001 mM as the product of the activity coefficient of [Ca²⁺]_ss and the concentration at the mouth of the channel. This assumes that once the channel is open, the concentration at the mouth does not change (Smith, 1996). This value was chosen because using the formulation presented by Luo-Rudy failed to give the initial fast rise and peak seen in experimentally measured L-type Ca²⁺ channels.

The final component of the L-type Ca²⁺ channel is the K⁺ current through the channel,

(38)

The variables y and {O + O_Ca} are the same as those described in Eq. 37, and the permeability P'_K is given in Eq. 1. The last term gives the driving force behind I_Ca,K.

The remaining equations for membrane currents and gating variables are taken from the Luo-Rudy Phase II model and are used to formulate differential equations (Appendix 1). It was necessary to rescale some of these currents: scaling factor of Na⁺-Ca²⁺ exchange k_NaCa = 5000 µA µF¹, maximum Na⁺-K⁺ exchange current I_NaK = 1.3 µA µF¹, maximum plateau K⁺ channel conductance _Kp = 0.00828 mS µF¹, maximum Na⁺ channel conductance _Na = 12.8 mS µF¹, maximum background Ca²⁺ current conductance _Ca,b = 0.006032 mS µF¹, and permeability of nonspecific current channel for K⁺ and Na⁺ P_ns(Ca) = 0.0 cms¹. The original Luo-Rudy Phase II values are shown in parentheses in Table 6 (Appendix 3). These mechanisms were rescaled as described above to maintain long-term homeostasis and improve AP shape.

	METHODS

Top Abstract Introduction Methods Results Concluding Remarks Appendix References

The full set of 30 ordinary differential equations was solved on a Silicon Graphics Indy workstation, using the Merson Modified Runge-Kutta 4th-Order Adaptive Step Algorithm (Kubicek and Marek, 1983), with a maximum step size of 0.1 ms and a maximum error tolerance of 10⁶. The error from all variables was normalized to ensure that all variables contributed equally to the global error calculation. These variables and their respective weights are V: 100 mV, [Na⁺]_i: 5 mM, [K⁺]_i: 140 mM, [Ca²⁺]_i: 0.001 mM, [Ca²⁺]_NSR: 2 mM, [Ca²⁺]_ss: 0.001 mM, and [Ca²⁺]_JSR: 20 mM. All other variables were given weights of 1.0. The results were visualized using Xmgr by Paul J. Turner and PV-WAVE by Precision Visuals.

The standard set of parameters and initial conditions used in the simulations is shown in the tables given in Appendix 3, unless specified in the figure legends or text. Because the original formulation for the Luo-Rudy model is for a ventricular cell of a small mammal like the guinea pig, when possible, the improved Ca²⁺ subsystem is also formulated using experimental data from small mammals. Action potentials were initiated by a 0.1 µA µF¹ current injection for 0.5 ms.

	RESULTS

Top Abstract Introduction Methods Results Concluding Remarks Appendix References

Single action potential

The model contains a new formulation for the L-type Ca²⁺ channel. The initial test of the model is to verify that the L-type Ca²⁺ current was consistent in shape and amplitude with experimental values measured by Grantham and Cannell (1996). The model was driven by the AP clamp used by Grantham and Cannell in their study of I_Ca (Fig. 5 A). Under control conditions and AP clamp, the model generates an I_Ca similar to experimentally measured values (Fig. 5 B, solid line). The contribution of Ca²⁺ inactivation is studied, following the experimental protocol (Grantham and Cannell, 1996) designed to deplete the SR. The external Ca²⁺([Ca²⁺]_o) is reduced to 0.1 mM, and the RyR open fraction (RyR_open) is set to 0.8 to mimic the effects of caffeine. The cell is then paced for 10 s at 1.0 Hz (free-running APs). The cell is returned to control conditions: no caffeine (RyR_open = P_O1 + P_O2) and [Ca²⁺]_o = 1.8 mM, and an AP clamp is applied. The I_Ca with depleted SR is shown by the dotted line (Fig. 5 B). The increase in I_Ca is due to the removal of Ca²⁺ inactivation. This is consistent with the findings of Grantham and Cannell (1996). The myoplasmic Ca²⁺ concentration ([Ca²⁺]_i) transient for the control under the AP clamp (Fig. 5 C; solid line) rises to 0.85 µM from a resting value of 0.10 µM. Under depleted SR conditions, the AP clamp yields a small [Ca²⁺]_i transient that rises to 0.12 µM from a resting value of 0.02 µM, indicating no Ca²⁺ release from the SR. The subspace Ca²⁺ ([Ca²⁺]_ss) in the control simulation under AP clamp (Fig. 5 D, solid line) rises to a peak of 22.2 µM. With a depleted SR, [Ca²⁺]_ss shows only a relatively small peak of 0.5 µM (Fig. 5 D, dotted line).

View larger version (38K): [in this window] [in a new window]   FIGURE 5   Simulations generated by the model in an AP clamp protocol (A-D) and in the free-running model (E-H). (A) Membrane potential trace from the AP clamp used by Grantham and Cannell (1996) to study ICa. (B) ICa generated by the model during the AP clamp protocol (control, solid line; depleted SR, dotted line; in all panels). The SR was depleted by pacing for 10 s at 1.0 Hz in low external Ca2+ (0.1 mM), with the RyR open fraction set to 0.8 to mimic the effects of caffeine. The external Ca2+ was returned to 1.8 mM and the first AP was recorded. (C) Myoplasmic Ca2+ transient. In the case with depleted SR (dotted line), the [Ca2+]i transient is attenuated. (D) The restricted subspace Ca2+ transient shows a higher peak and faster dynamics than [Ca2+]i. (E) AP generated by the free-running model. The AP with a depleted SR is longer than the control AP. (F) ICa in free-running model. (G) Myoplasmic Ca2+ transient. (H) Subspace Ca2+ transient. Parameters and initial conditions are given in the tables in Appendix 3.

The protocol described above is applied to the free-running model to yield a typical action potential. The initial conditions for this simulation are obtained during diastole from a simulation of a cardiac myocyte paced at 1.0-Hz stimulus frequency. The AP is triggered at 10 ms. The control AP rises to a peak value of 50.4 mV (Fig. 5 E, solid line). The control AP is shorter in duration and has a lower plateau than the AP with depleted SR (Fig. 5 E, dotted line). The I_Ca generated by the model (Fig. 5 F, solid line) is similar to those observed in experiments (Grantham and Cannell, 1996) and very close to those observed by Puglisi and Bers (personal communication). The L-type Ca²⁺ current quickly activates, reaching its peak of 6.46 µA µF¹ within 2 ms of the stimulus. The initial decrease in the current, giving rise to its spike-like appearance, results from a decline in the L-type Ca²⁺ current as V rises to the peak of the AP. As Fig. 5 B shows, at 50 mV the peak L-type Ca²⁺ current through the channel is greatly diminished. A secondary reduction of the current after ~30 ms is produced by Ca²⁺ inactivation. Later in the AP, I_Ca is further diminished because of voltage-dependent inactivation (~100 ms). At this time, the high value of [Ca²⁺]_ss drives the channel to mode Ca, which has infrequent openings and hence lowers open probability. With SR depletion, I_Ca increases slightly in amplitude and duration (Fig. 5 F, dotted line) because of the removal of Ca²⁺ inactivation. The longer duration and higher plateau of the SR-depleted AP is due to these changes in I_Ca. The [Ca²⁺]_i transients rises to 0.84 µM from a resting value of 0.10 µM (Fig. 5 G, solid line). Under depleted SR conditions, the AP yields a small [Ca²⁺]_i transient that rises to 0.12 µM from a resting value of 0.02 µM, indicating no Ca²⁺ release from the SR (Fig. 5 G, dotted line).

The Ca²⁺ transient in the control simulation (Fig. 5 G, solid line) is similar in magnitude and time course to experimental data from the guinea pig (Isenberg, 1995). A typical Ca²⁺ transient measured experimentally reaches a peak of 1.0 µM from a resting value of 0.1 µM. In the model at 1.0 Hz pacing frequency, [Ca²⁺]_i rises from 0.10 µM during diastole to a peak value of 0.84 µM during an AP, consistent with these estimates of changes in bulk [Ca²⁺]_i. In the guinea pig ventricular myocyte, this Ca²⁺ transient has a duration of ~0.3 s and has a tonic elevation that follows the peak systolic [Ca²⁺]_i (Isenberg, 1995). The model reproduces this tonic elevation by the shoulder seen in Fig. 5 G to the right of the peak. In experiments, Ca²⁺ reaches its peak value (latency) within 20-40 ms and decays within 80-120 ms (Isenberg, 1995). In the model, the latency will not exceed this limit through a wide variety of pulse protocols. For example, in Fig. 5 G the latency is 10.0 ms. In contrast, the Ca²⁺ in restricted subspace ([Ca²⁺]_ss) cannot yet be quantified experimentally, but is thought to be at least an order of magnitude larger than bulk systolic [Ca²⁺]_i (Santana et al., 1996). In the model, [Ca²⁺]_ss reaches a peak value of 22.3 µM, as shown in Fig. 5 H, starting from a resting value of 0.14 µM. Note that [Ca²⁺]_ss shows a different time course, with much steeper onset and decay, than [Ca²⁺]_i. This demonstrates that [Ca²⁺]_ss and [Ca²⁺]_i are not just scaled versions of each other, because there is a very different set of sources and sinks for Ca²⁺ in each of the compartments. These will now be examined in more detail.

The intracellular Ca²⁺ fluxes cannot currently be measured, so the following simulation results cannot be compared to experimental data. However, these model results reveal dynamics of the underlying Ca²⁺ subsystem, and predict how these signals might appear in real cells. The RyR flux (J_rel) has a peak value ~40 times larger than that of the J_ICa, consistent with the large amplification estimated for CICR (Cannell et al., 1994). When total Ca²⁺ influx is summed over the course of the AP, ~10 times more Ca²⁺ enters the subspace from the RyRs than from the L-type channels, consistent with a number of estimates (Isenberg, 1995; Bers, 1991).

Efflux from the subspace is by diffusion to the myoplasm, which is described by the transfer flux out of the subspace (J_xfer). The transfer time constant _xfer is set to a value that is a compromise between two extremes. Small values produce a large concentration peak in the subspace, slow transfer to bulk myoplasm, and possibly restimulate CICR. Large transfer rates cause the influx of trigger Ca²⁺ through the L-type channel to dissipate quickly. This leads to a long latency for the Ca²⁺ transient due to the slowed initiation of CICR. In extreme cases, the rise in concentration of the restricted subspace may be ineffective in producing regenerative CICR.

The major source of influx to the bulk myoplasmic space is from J_xfer, the transfer from the subspace to the bulk myoplasm. Another small source of flux into the myoplasm is from the leak of Ca²⁺ out of the SR. The major efflux is the uptake into the SR by the Ca²⁺-ATPase (J_up). A much smaller source of efflux is the Na⁺/Ca²⁺ exchanger (I_NaCa), but only very late in the AP; this exchanger works in reverse earlier in the AP, producing an influx of Ca²⁺. I_NaCa is electrogenic; in forward mode, three Na⁺ ions enter for every Ca²⁺ ion extruded, producing a net inward current. Because of this electrogenic nature, the exchanger shows voltage dependence and works in reverse mode at depolarized potentials. The switch from Na⁺ extrusion to Ca²⁺ extrusion occurs later in the AP than in the original Luo-Rudy formulation; however, this late switch from reverse to forward mode is consistent with the I_NaCa measured by Grantham and Cannell (1996). The reason for the difference is due to the altered shape of the [Ca²⁺]_i transient. The Luo-Rudy phase II model has a triangular calcium transient, whereas the Ca²⁺ transients generated by our model more closely simulate experimentally determined Ca²⁺ transients. The smallest source of efflux is the sarcolemmal Ca²⁺ pump. This pump has a high affinity and low throughput, so that it generally contributes little efflux under normal conditions. Its contribution may be important, however, in long-term Ca²⁺ homeostasis or during pathological conditions.

All uptake of Ca²⁺ (J_up) is into the network SR (NSR) compartment (Fig. 7 B, dashed line). There is an initial rise in [Ca²⁺]_NSR resulting from the increase in J_up when [Ca²⁺]_i rises during the AP. All SR release (J_rel) comes from the junctional SR (JSR). This release produces the drop in [Ca²⁺]_JSR. The depletion of this store is one factor in terminating CICR. Another factor terminating CICR is the decrease in open fraction of the RyRs as a result of adaptation. The JSR store is replenished from a transfer flux from NSR (J_tr).

To understand the relative contributions of RyR adaptation and JSR depletion to termination of Ca²⁺ release, we performed simulations that varied either the rate of adaptation or level of JSR Ca²⁺ depletion (Figs. 6 and 7). The first set of simulations studied the contribution of adaption to the termination of Ca²⁺ release from the SR. To eliminate adaptation without perturbing the system too drastically, we set the rate constant k_c⁺ for migration into the RyR adapted state C₂ to zero (Fig. 4) and paced the cell at 1.0 Hz until a stable set of APs was reached. When compared to the control simulation (k_c⁺ = 0.018 ms¹; solid line), this results in almost no decrease in AP duration (APD), almost no change in Ca²⁺ release (Fig. 6 A, dotted line), and a slightly slower time course for [Ca²⁺]_JSR refilling (Fig. 6 B, dotted line). Fig. 6 C (dotted line) shows that whereas the opening of the RyR is not significantly delayed, the closing is delayed slightly, resulting in a small increase in the time interval during which RyR open probability is high. The small magnitude of the changes in Fig. 6 show that for the most part, model responses are similar before and after elimination of RyR adaptation.

View larger version (21K): [in this window] [in a new window]   FIGURE 6   The role of RyR adaptation in the termination of Ca2+ release during an AP: 1) control simulation with the standard parameter set (kc+ = 0.018 ms1; solid line); 2) no adaptation (kc+ = 0.0 s1; dotted line); and 3) fast adaptation (kc+ = 0.072 ms1; dashed line). The AP with no adaptation is almost indistinguishable from the control. With fast adaptation, the AP has a slightly higher plateau and slightly longer duration. (A) The [Ca2+]i transient with no adaptation is similar to control, whereas with fast adaptation the latency is increased and the amplitude is reduced. (B) Depletion of [Ca2+]JSR is delayed and reduced in the case of fast adaptation. In the case of no adaptation, [Ca2+]JSR depletion recovers slightly more slowly than control. (C) With no adaptation, the RyR open fraction transient lasts slightly longer than the control. The RyRopen transient is smaller, is of shorter duration, and has delayed onset compared to control.

View larger version (22K): [in this window] [in a new window]   FIGURE 7   The role of SR depletion on the termination of Ca2+ release during an AP: 1) control simulation (solid line); 2) SR clamp ([Ca2+]JSR = 1.175 mM) at resting value ([Ca2+]JSR = 1.175 mM; dotted line); 3) partial SR clamp at 50% resting value ([Ca2+]JSR = 0.588 mM; long dashed line). The dashed line shows [Ca2+]NSR in the control simulation. The action potential in the case of SR clamp does not converge, because of never-ending release from the SR (A). The partial SR clamp yields an AP with elevated plateau and longer duration. (A) The myoplasmic Ca2+ transient with SR clamp does not terminate. With the partial SR clamp, there is no Ca2+ release from the SR. (B) [Ca2+]JSR varies from a resting value of 1.175 mM to a minimum value of 0.108 mM during a control AP. With the SR clamp the value remains fixed at 1.175 mM, whereas during the partial SR clamp the value remains constant at 0.588 mM. (C) During SR clamp conditions the RyR open fraction remains at values similar to the peak values during a control. During partial SR clamp, RyRopen barely rises at all, indicating minimal release from the SR.

In the next set of simulations, we examined the effects of increasing the rate of RyR adaptation on Ca²⁺ release. The parameter k_c⁺ was increased from 0.018 ms¹ (Fig. 6, solid line) to 0.072 ms¹ (Fig. 6, dotted line), effectively increasing the rate of adaptation by 4.0 times. Faster rates of adaptation were not used because they resulted in alternans. The model was paced at 1.0 Hz until a stable set of APs was obtained. With the faster adaptation rate, the AP is only slightly longer than the control, with a very slightly elevated plateau (not shown). This lengthening is due to an increase in I_Ca, caused by a reduction of Ca²⁺ inactivation in response to the smaller Ca²⁺ transients in the subspace, which results in a small Ca²⁺ transient in the myoplasm. In addition to the smaller magnitude of the [Ca²⁺]_i transient, the transient also shows a longer latency to peak than control (solid line). The longer latency is a result of a lower gain for CICR. This occurs because more of the RyRs are in the adapted state (P_C2), and thus there is a smaller fraction that can quickly make a transition from P_C1 to P_O1 in response to Ca²⁺ and hence contribute to regenerative release. Depletion of [Ca²⁺]_JSR has a delayed onset and a smaller magnitude than in the control (Fig. 6 B). With fast adaptation, the resting [Ca²⁺]_JSR is a little lower than the control levels. The open fraction transient has decreased magnitude and duration (Fig. 6 C) compared to control. Hence the simulations suggest that RyR adaptation has a small effect on modulating the Ca²⁺ transient in a single AP.

The next set of simulations shows the effects of SR load on SR release termination. Here [Ca²⁺]_JSR is clamped by reducing its time derivative to zero (SR clamp). With this SR clamp (Fig. 7, dotted line), the system becomes unstable and Ca²⁺ release from the SR does not terminate. This is evident in the AP increasing without bound (not shown) due to increased Na⁺-Ca²⁺ exchange, which results from a greatly increased [Ca²⁺]_i (Fig. 7 A). The level of [Ca²⁺]_JSR remains constant at its clamped value of 1.175 mM (Fig. 7 B). RyR_open for the SR clamp (Fig. 7 C, dotted line) rises to the peak value seen in the control simulation (solid line) and remains there. Even in the case of fast adaptation described above, Ca²⁺ release does not terminate under the conditions of SR clamp (not shown). These simulations suggest that JSR depletion is essential for termination of Ca²⁺ release from the SR. Upon first glance, it might seem that in these simulations Ca²⁺ release termination by the SR depletion is due to a lack of releasable Ca²⁺ in the JSR. This, however, is not the case. Instead, it is the reduction of the gradient for Ca²⁺ between the JSR and the subspace that terminates release. Once this gradient diminishes, there is not enough Ca²⁺ entering the subspace to trigger and sustain regenerative release. To demonstrate this point, [Ca²⁺]_JSR is clamped to half its resting value (0.588 mM; Fig. 7, long dashed line). With this partial SR clamp, the AP had longer duration with an elevated plateau compared with the control (not shown). There is no regenerative release from the SR, as indicated by the small [Ca²⁺]_i transient (Fig. 7 A) and the tiny RyR open fraction (Fig. 7 C). Even a large (10 µM) step in [Ca²⁺]_ss at stimulation does not initiate CICR. However, if P_O1 is increased transiently to 0.7 (and P_C1 decreased accordingly) at stimulation, CICR does occur, depleting the [Ca²⁺]_JSR to levels similar to those seen during the control AP. Thus the Ca²⁺ release flux from the SR is insufficient to maintain CICR, even though there is still releasable Ca²⁺ in the JSR. During a control AP, the level of [Ca²⁺]_JSR falls below the value of the partial SR clamp (Fig. 7 B), indicating that there is still releasable Ca²⁺ in the store. From these sets of simulations, it is clear that JSR Ca²⁺ depletion terminates CICR, whereas RyR adaptation has only a small effect on Ca²⁺ release. Thus, during the course of an AP, Ca²⁺ release from the SR is terminated by depletion of the JSR, whereas RyR adaptation only modulates Ca²⁺ release. Furthermore, the amount of Ca²⁺ released via RyR during an AP can be assessed to show that the SR as a whole is not depleting most of its Ca²⁺ content. In the model the SR is subdivided into two functional compartments, the JSR (Fig. 7 B, solid line) and NSR (Fig. 7 B, dashed line). Calculations show that the amount of Ca²⁺ that leaves the JSR is ~28% of the total SR Ca²⁺ content with 1.8 mM external Ca²⁺. This compares well with the value of 35% by Bassani and co-workers (Bassani et al., 1995) for an external Ca²⁺ concentration of 2.0 mM measured in ferret ventricular myocytes. According to the trend they observed, a lower external Ca²⁺ will yield a lower fraction released.

Pacing protocols

The results above illustrate model behavior over the relatively short time scale of one AP. However, many phenomena associated with E-C coupling depend on the stimulus interval, and hence evolve over a time course of many APs. These phenomena, known as the interval-force relations, describe the changes in force generation as pacing frequency changes. In most species, when the pacing frequency is increased, the amount of force generated in response to an AP rises to a new steady value after a transient decrease (Bers, 1991). Upon a reduction of pacing frequency to its original value, the force generated returns to its previous steady level after a transient increase. Sample results for the step changes in pacing frequencies generated by the model are shown in Fig. 8. In this protocol, the stimulus frequency is transiently increased from 0.5 Hz to 1.5 Hz and then returned to 0.5 Hz. The pacing at 0.5 Hz is initially maintained for a long period, so that the simulated myocyte is producing a stable output, as shown by the first seven beats in Fig. 8. The amplitude of the AP decreases slightly with the increase in stimulation rate (Fig. 8 A). The effects of different pacing frequencies upon APs will be described later (Figs. 8-12). When the stimulus frequency is raised to 1.5 Hz, the Ca²⁺ transient decreases initially, but then rises to a new plateau level over the course of the next few beats (Fig.