HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 286/3/H1154    most recent 00168.2003v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (12) Citing Articles via Google Scholar Google Scholar Articles by Bondarenko, V. E. Articles by Rasmusson, R. L. Search for Related Content PubMed PubMed Citation Articles by Bondarenko, V. E. Articles by Rasmusson, R. L.

A model of graded calcium release and L-type Ca²⁺ channel inactivation in cardiac muscle

Department of Physiology and Biophysics, University at Buffalo, State University of New York, Buffalo, New York 14214

Submitted 20 February 2003 ; accepted in final form 6 November 2003

	ABSTRACT

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX A: DIFFERENTIAL... APPENDIX B: CALCIUM CURRENTS APPENDIX C: CALCIUM FLUXES REFERENCES

 
We have developed a model of Ca²⁺ handling in ferret ventricular myocytes. This model includes a novel L-type Ca²⁺ channel, detailed intracellular Ca²⁺ movements, and graded Ca²⁺-induced Ca²⁺ release (CICR). The model successfully reproduces data from voltage-clamp experiments, including voltage- and time-dependent changes in intracellular Ca²⁺ concentration ([Ca²⁺]_i), L-type Ca²⁺ channel current (I_CaL) inactivation and recovery kinetics, and Ca²⁺ sparks. The development of graded CICR is critically dependent on spatial heterogeneity and the physical arrangement of calcium channels in opposition to ryanodine-sensitive release channels. The model contains spatially distinct subsystems representing the subsarcolemmal regions where the junctional sarcoplasmic reticulum (SR) abuts the T-tubular membrane and where the L-type Ca²⁺ channels and SR ryanodine receptors (RyRs) are localized. There are eight different types of subsystems in our model, with between one and eight L-type Ca²⁺ channels distributed binomially. This model exhibits graded CICR and provides a quantitative description of Ca²⁺ dynamics not requiring Monte-Carlo simulations. Activation of RyRs and release of Ca²⁺ from the SR depend critically on Ca²⁺ entry through L-type Ca²⁺ channels. In turn, Ca²⁺ channel inactivation is critically dependent on the release of stored intracellular Ca²⁺. Inactivation of I_CaL depends on both transmembrane voltage and local [Ca²⁺]_i near the channel, which results in distinctive inactivation properties. The molecular mechanisms underlying many I_CaL gating properties are unclear, but [Ca²⁺]_i dynamics clearly play a fundamental role.

myocytes; computer simulation; inactivation; calcium-induced calcium release; excitation-contraction coupling

Recently, CICR has been observed at the subcellular level in the form of Ca²⁺ "sparks," which are extremely large, extremely localized increases in [Ca²⁺] (18, 19, 29, 34, 35). Spontaneously occurring sparks have been observed. They can trigger a subcellular Ca²⁺ wave that can propagate throughout the cell (17, 30, 35). Ca²⁺ sparks were initially thought to reflect the opening and closing of a single SR ryanodine receptor (RyR) (18). However, they have been subsequently shown to be the result of the opening of a single L-type Ca²⁺ channel (38), which is sufficient to initiate CICR from a cluster of RyRs (36). The number and timing of sparks depend in a nonmonotonic manner on the amplitude of the membrane depolarization.

One of the most important features of intracellular [Ca²⁺] cycling is the graded release of Ca²⁺ from the SR (9, 11, 14). Recently, Rice et al. (57) published a model of graded CICR based on a Monte-Carlo simulation with a bell-shaped profile for the both integrated and peak RyR Ca²⁺ release fluxes. Our model generates graded release through biophysically distinct mechanisms. The Rice et al. model (57) postulates that graded release occurs through stochastic changes in the release channel. Subsequent variants (64) added additional deterministic elements of structure, but the main mechanism of graded release remained stochastic. In contrast, our model assumes a physiological variability in diadic junction morphology to generate graded release. The models are, therefore, inherently different at a qualitative level, and there will be testable differences in predicted behavior between these models.

In this paper, we present a new model of cardiac Ca²⁺ handling that incorporates a novel Markov model for L-type Ca²⁺ channels and graded SR Ca²⁺ release in response to voltage-clamp simulations. We hypothesize that the structural inhomogeneity of the local calcium release subsystems is responsible for graded release. To model this situation, we introduced eight calcium release subsystems with different characteristics. The model myocyte has a series of physically distinct subsystems that represent the area where the junctional SR comes in close proximity to the T-tubular membrane. Each subsystem is a relatively small subspace volume bounded by the sarcolemma, which contains between one and eight L-type Ca²⁺ channels, and the junctional SR. The SR surface exposed to the subsystem has a cluster of RyRs that release Ca²⁺ into the subsystem volume. The subsystems are distributed throughout the model cell, with the number of L-type Ca²⁺ channels in each subsystem determined by a binomial weighting function. CICR is triggered from the cluster of RyRs in the subsystem when the amount of Ca²⁺ entering the subspace via the L-type Ca²⁺ channels is sufficiently large. We also developed a new Markov model for the L-type Ca²⁺ channel. This L-type Ca²⁺ channel current (I_CaL) gating model is based on emerging structure and function data on L-type channels and related voltage-gated Na⁺ and K⁺ channels. It contains one activation pathway with four closed states and a final open state. Our model explicitly includes Ca²⁺- and voltage-dependent inactivation mechanisms and recovery pathways from each inactivated state. The model reproduces the experimental behavior of I_CaL in isolated ventricular cardiac myocytes. It also offers an explanation of relationship between membrane voltage and the number of elementary Ca²⁺ sparks observed (as determined by confocal scanning microscopy experiments).

	METHODS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX A: DIFFERENTIAL... APPENDIX B: CALCIUM CURRENTS APPENDIX C: CALCIUM FLUXES REFERENCES

 
Model myocyte. The currents, fluxes, and physical layout of the model cell are shown in Fig. 1. Depolarization of the cell membrane opens voltage-gated L-type Ca²⁺ channels, which allow Ca²⁺ to enter the subspace volume (V_ss) and increase the subspace Ca²⁺ concentration ([Ca²⁺]_ss). The entering calcium (I_CaL) triggers CICR from the junctional SR (J_rel), which induces further CICR and inactivates I_CaL (65). Ca²⁺ in the subspace diffuses from the subsystem to the general cytosol (J_xfer), where it contributes to [Ca²⁺]_i. Calcium leaves the general cytosol by being pumped back into the SR by SR Ca²⁺-ATPase (J_up), pumped across the cell membrane by the sarcolemmal Ca²⁺ pump [Ca²⁺ pump current (I_pCa)], or extruded from the cell by the Na⁺/Ca²⁺ exchange mechanism [Na⁺/Ca²⁺ exchange current (I_NaCa)]. Calcium buffers present in the general cytosol also remove Ca²⁺ from free cytosolic activity: Ca²⁺ is modeled as being at instantaneous equilibrium with calmodulin, but there are fast and slow time-dependent binding sites for troponin (J_trpn). Calcium that is pumped into the network SR diffuses to the junctional SR (J_tr), which can store large quantities of Ca²⁺ due to the presence of calsequestrin. Calcium is assumed to be in instantaneous equilibrium with calsequestrin. A small leak current takes Ca²⁺ directly from the network SR to the general cytosol (J_leak), and the background Ca²⁺ current (I_Cab) transfers Ca²⁺ from the extracellular space to the main cytosol.

View larger version (31K): [in this window] [in a new window]   Fig. 1. A: spatial organization and Ca2+ fluxes in the ferret ventricular myocyte model. Ca2+ enters the cell via L-type Ca2+ channel current (ICaL) and initiates Ca2+-induced Ca2+ release (CICR) from the junctional sarcoplasmic reticulum (JSR), Jrel. Ca2+ diffuses from the subcellular compartment to the main cytosol, Jxfer, where it binds to calmodulin and troponin, Jtrpn. Ca2+ is removed from the cytosol by sarcolemmal Na+/Ca2+ exchanger current (INaCa) and Ca2+ pump current (IpCa) or is pumped into the SR via SR Ca2+-ATPase, Jup. There is a Ca2+ background current (ICab) into the cytosol. Ca2+ in the network SR (NSR), [Ca2+]NSR, mainly diffuses to the JSR, Jtr, but a small fraction leaks back to the cytosol, Jleak. Ca2+ in the JSR, [Ca2+]JSR, binds to calsequestrin. [Ca2+]o, extracellular Ca2+ concentration; [Ca2+]i, intracellular Ca2+ concentration. B: schematic drawing showing a cross section of 3 of the 8 possible Ca2+ release subsystems, with 5, 1, or 8 L-type Ca2+ channels. Ca2+ entering through the L-type Ca2+ channels initiates CICR from the group of ryanodine receptors (RyRs) on the opposing SR, resulting in a Ca2+ flux, Jrel,k, from the JSR volume (VJSR) into the subspace volume (Vss,k). Ca2+ diffuses from the subspace volume to the cytosol, producing the flux Jxfer,k. DHPR, dihydropyridine receptor.

Differential equations. The model consists of 117 time-dependent differential equations (APPENDICES A–C), which are solved using a fourth-order Runge-Kutta method with a time step of 0.0001 ms. This was programmed in Fortran 90 and run on a 500-MHz DEC Alpha Workstation under Digital UNIX version 4. For simulation of a single voltage-clamp trace, 1 s of data required 400 s of computational time. The initial conditions for all variables are detailed in Table 1. The model is designed to simulate voltage-clamp experiments and so has only four transmembrane currents: I_CaL, I_Cab, I_pCa, and I_NaCa. Neither intracellular Na⁺ concentration ([Na⁺]_i) nor intracellular K⁺ concentration are modeled, and there are no transmembrane Na⁺ or K⁺ currents

(1)

where I_total is total current.

Calcium release subsystem and SR calcium flux. The model has physically distinct subsystems that represent the small intracellular spaces into which transmembrane Ca²⁺ flux through the L-type Ca²⁺ channel enters, CICR is initiated, and SR Ca²⁺ is released. Each subsystem has between one and eight L-type Ca²⁺ channels and a cluster of RyRs. Figure 1B shows three separate Ca²⁺ release subsystems, Ca_ss5, Ca_ss1, and Ca_ss8, with five, one, and eight L-type Ca²⁺ channels, respectively, and an I_CaL magnitude given by I_CaL,5, I_CaL,1, and I_CaL,8, respectively. The model has a binomial distribution of subsystems (Ca_ss1 to Ca_ss8) as a default. However, we also tested the model with three other weighting functions (uniform, exponential, and Maxwellian) to determine whether the particular weighting function used affected graded release. All distributions were normalized to produce a mean value of four L-type Ca²⁺ channels per RyR cluster. The weighting functions for L-type Ca²⁺ channel distributions (bin_k, unif_k, exp_k, and maxw_k) are described in Table 2.

The choice of eight subsystems was based on estimates of RyR cluster size and the ratio of dihydropyridine receptors (DHPR) to RyRs (8, 28). We assume that the average cluster contains 40 RyRs and 4 DHPRs, which gives a ratio of 10:1 (8). A single RyR cluster covers a surface area of 360 nm x 360 nm and is separated from the myocyte membrane by a gap of 12 nm (28). A whole cell current of 1 nA is the result of 5,000 Ca²⁺ channels opening (6), so the average number of active release subsystems is 1,250, with 4 DHPRs per subsystem. The total subspace volume V_ss is therefore estimated to be 2.0 x 10^–9 µl (360 nm x 360 nm x 12 nm x 1,250). There is a linear relationship between V_ss,k and the number of L-type Ca²⁺ channels (DHPRs) in our model. The model default is to have 40 RyRs per release subsystem, independent of the size of the subsystem. However, we also tested the model with a variable number of RyRs per subsystem where the number of RyRs was a factor of 10 greater than the number of DHPRs (i.e., 10–80 RyRs per subsystem with 1–8 L-type Ca²⁺ channels).

For the subsystems with k L-type Ca²⁺ channels, the L-type Ca²⁺ current is I_CaL,k, the Ca²⁺ concentration in V_ss,k is [Ca²⁺]_ss,k, the CICR flux is J_rel,k, and the diffusion of Ca²⁺ from the subspace volume to the main cytosol is J_xfer,k. V_ss,k depends linearly on the number of L-type Ca²⁺ channels present. The time constant of the transfer of Ca²⁺ from the subspace to the main cytoplasm (_xfer,k) depends linearly on size of the subsystem volume (see Table 2). The junctional SR volume (V_JSR) is constant for all subsystems (see Table 3), although the Ca²⁺ efflux from the SR, J_rel,k, is different for each subsystem.

The total Ca²⁺ flux depends on the sum of the activity in all subsystems with k L-type Ca²⁺ channels and is described by the following equation

(2)

where x is rel for release, tr for transfer from the junctional SR to the network SR, and xfer for transfer from the subspace to the cytoplasm.

The SR release flux (J_rel,k) depends on the difference in [Ca²⁺] between the junctional SR ([Ca²⁺]_JSR,k) and the subspace volume ([Ca²⁺]_ss,k), the sum of probabilities for opening of the RyRs (P_O1,k and P_O2,k), and the maximum permeability of RyR channels (v₁)

(3)

The magnitude of the Ca²⁺ flux from the network to junctional SR (J_tr,k) depends on the differences in [Ca²⁺] (i.e., [Ca²⁺]_NSR and [Ca²⁺]_JSR,k) and the time constant of diffusion (_tr)

(4)

Ca²⁺ flux from the subspace volume to the cytoplasm (J_xfer,k) depends on the relative concentrations of calcium in the subspace ([Ca²⁺]_ss,k) and the cytosol ([Ca²⁺]_i) and the time constant of diffusion (_xfer,k) from the subspace volume (V_ss,k) to the general cytosol (V_myo)

(5)

The uptake Ca²⁺ flux to the SR (J_up) depends on [Ca²⁺]_i, the maximum SR Ca²⁺-ATPase pump rate (v₃), and the half-saturation constant for SR Ca²⁺-ATPase (K_m,up)

(6)

The rate at which calcium binds to troponin (J_trpn) is

(7)

where [LTRPN]_tot is the total concentration of the myoplasmic troponin low-affinity site, [HTRPN]_tot is the total concentration of the myoplasmic troponin high-affinity site, [LTRPNCa] is the concentration of low-affinity troponin-binding sites with Ca²⁺ bound, [HTRPNCa] is the concentration of high-affinity troponin-binding sites with Ca²⁺ bound, is the Ca²⁺ on rate for troponin high-affinity sites, is the Ca²⁺ off rate for troponin high-affinity sites, is the Ca²⁺ on rate for low-affinity troponin sites, and is the Ca²⁺ off rate for low-affinity troponin sites (Table 4). Calmodulin, the other cytosolic Ca²⁺ buffer, is assumed to be in instantaneous equilibrium with [Ca²⁺] (see Eqs. 15, 16, 19, and 20).

The leak of Ca²⁺ from the SR to the cytosol (J_leak) depends on the relative concentrations of Ca²⁺ and the leak rate constant (v₂)

(8)

L-type Ca²⁺ channels. Modeling the inactivation behavior of the L-type Ca²⁺ channel has long been a difficulty. Its complex kinetics involve both voltage and calcium dependence. Reproducing this behavior is a current research topic with different competing models (31, 39, 64, 68). One striking piece of information to emerge from the field of molecular biology is the sequence homology between Ca²⁺ channels and Na⁺ or K⁺ channels. Several model structures can reproduce various aspects of native channel gating (31, 48, 68). Our choice was to assume that gating of Ca²⁺ channels was similar to gating of Na⁺ and K⁺ channels. In these channels, activation is rapid and intrinsically voltage sensitive. In contrast, inactivation has multiple mechanisms, all of which are intrinsically voltage insensitive and coupled to activation (for a review, see Ref. 56). In our case, there are two obvious mechanisms of inactivation, fast Ca²⁺-dependent inactivation and slow Ca²⁺-insensitive inactivation. Fast calcium- and calmodulin-sensitive inactivation may be similar to N-type inactivation (48) and slow Ca²⁺-independent inactivation may be similar to C-type inactivation (77). This formulation can reproduce the behavior of the Ca²⁺ channel but also predicts that mutagenesis experimental results will parallel those for Na⁺ and K⁺ channels. Furthermore, models of open channel and use-dependent blockers are strongly sensitive to the exact formulation used to describe ion channel gating (37) so that predicted behavior of use-dependent drugs may be strongly influenced by this model formalism.

A new Markov model was developed for the L-type Ca²⁺ channel with four closed states (C₁, C₂, C₃, and C₄), one open state (O), and three inactivated states (I₁, I₂, and I₃). The states and transitions are shown in Fig. 2. Inactivation in this model is based on structural and functional homology with voltage-gated Na⁺ and K⁺ channels. The core "voltage-dependent" inactivation that persists in the absence of Ca²⁺-dependent inactivation is thought to be similar to K⁺ channel C-type inactivation (77) or the slow inactivation in Na⁺ channels (2, 72). The O I₂ transition represents the channel undergoing "voltage-dependent" inactivation. Fast Ca²⁺-dependent inactivation arises from a tethered calmodulin-binding ancillary subunit (47, 48). This subunit is modeled as a "ball-and-chain" type of inactivation (76) in which the ball requires Ca²⁺ to bind to "activate" it (20). The O I₁ transition corresponds to Ca²⁺-dependent inactivation, which is controlled by the Ca²⁺-dependent rate constant _k. I₃ represents a channel that is both voltage and Ca²⁺ inactivated. In a similar fashion to Na⁺ and K⁺ channels, activation is required before inactivation can occur (for a review, see Ref. 26). Direct transitions between the C₄ closed state and the three inactivated states were introduced into the model to reproduce voltage-dependent recovery from inactivation (27). All rate constants were chosen to satisfy thermodynamic equilibrium, i.e., the product of the rate constants in any loop of channel states in a clockwise direction is equal to the counterclockwise product.

View larger version (16K): [in this window] [in a new window]   Fig. 2. State diagram for the L-type Ca2+ channel. There are four closed states (C1–C4), one open state (O), and three inactivated states (I1–I3). Transitions between the closed states and the open state are voltage dependent. Transition to the I1 inactivated state is Ca2+ dependent and that to I2 is voltage dependent. The third inactivated state, I3, represents a state that is both voltage- and Ca2+-dependent inactivation. Transitions from the inactivated states I1--I3 to the closed state C4 correspond to recovery from inactivation. Kpcf and Kpcb, voltage-insensitive rate constants; k, Ca2+-dependent rate constant.

Each Ca²⁺ subsystem contains a different number of L-type Ca²⁺ channels and different subspace Ca²⁺ concentrations, [Ca²⁺]_ss,k. The calcium-induced inactivation of the channel, and therefore the behavior of the L-type Ca²⁺ channel, will be different in each subsystem. The total I_CaL is the sum of all the subsystem Ca²⁺ currents

(9)

where N is the number of calcium release subsystems and I_CaL,k is I_CaL in a subsystem with k channels, given by

(10)

where is the specific maximum conductivity of I_CaL, V is membrane potential, E_Ca,L is the channel reversal potential, and O_k is the probability of the channel being open.

Other Ca²⁺ currents. The model includes three other Ca²⁺ transmembrane pathways: I_pCa, I_NaCa, and I_Cab. I_pCa is based on that of Luo and Rudy (39)

(11)

where is the maximum I_pCa and K_m,pCa is the Ca²⁺ half-saturation constant. I_Cab is modeled as a simple linear ohmic conductor

(12)

(13)

where G_Cab is background Ca²⁺ conductance, E_CaN is the reversal potential, F is Faraday's constant, T is the absolute temperature (modeled at 295 K, i.e., room temperature), and R is the ideal gas constant.

I_NaCa has a stoichiometry of 3 Na⁺:1 Ca²⁺, so there is a net inward movement of a single charge associated with the efflux of a single calcium ion. The voltage dependence of the exchanger was calculated as the movement of a charge through a single peak energy barrier. The equation representing I_NaCa is based on those of Refs. 42 and 43, which were adapted by Luo and Rudy (39)

(14)

where k_NaCa is the exchange mechanism scaling factor, K_m,Na is the Na⁺ half-saturation constant, k_m,Ca is the Ca²⁺ half-saturation constant, K_sat is the Na⁺/Ca²⁺ exchange saturation factor at very negative potentials, and is the partition parameter, which indicates the relative position of the energy minima in the membrane (Table 5). Na⁺ dynamics are not included in the model, so [Na⁺]_i and extracellular Na⁺ concentration ([Na⁺]_o) are constants, which are given in Table 6.

Calcium concentrations. The equations governing the rate of change of [Ca²⁺] in the various cellular compartments (e.g., V_ss,k, V_JSR, V_NSR, and V_myo) are as follows

(15)

(16)

(17)

(18)

where

(19)

(20)

(21)

t is time, A_cap is the capacitative membrane area, V_myo is the myoplasmic volume, V_NSR is the network SR volume, [CMDN]_tot is the total concentration of myoplasmic calmodulin, is the Ca²⁺ half-saturation constant for calmodulin, [CSQN]_tot is the network SR calsequestrin concentration, and is the Ca²⁺ half-saturation constant for calsequestrin. Note that the fluxes and variables with a k index in Eqs. 15–17 refer to the kth Ca²⁺ release subsystem. The weighted sum of the subsystem contributions are denoted by the same symbol but lack a k index (see Eq. C1 in APPENDIX C).

The rate of change of occupancy of low- and high-affinity troponin-binding sites ([LTRPNCa] and [HTRPNCa], respectively) is given by the equations

(22)

(23)

Ca²⁺ binding to calsequestrin and calmodulin is considered to be instantaneous (31). In contrast, Ca²⁺ binding to troponin is much slower and is therefore described by time-dependent differential equations (31). The validity of the rapid buffering approximation (70) under physiological conditions has been verified by Smith et al. (63). Smith et al. (63) showed that approximating calmodulin buffering as an instantaneous event is appropriate for calmodulin but not for troponin. The binding kinetics of troponin can deviate from an instantaneous approximation by an order of magnitude at a time scale of 1 ms.

Each of the eight calcium subsystems can have a different value for [Ca²⁺]_ss,k, so the state of the L-type Ca²⁺ channels will vary between subsystems. The time course of the channel for the conducting state (O_k), the closed states (C_1,k–C_4,k), and the inactivated states (I_1,k–I_3,k) for the kth Ca²⁺ release subsystem is determined by the following equations

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

where K_pcf and K_pcb are constants (listed in Table 7) and

(32)

(33)

(34)

where K_pc,max is the maximum time constant for Ca²⁺-induced inactivation and K_pc,half is the half-saturation constant for Ca²⁺-induced inactivation (listed in Table 7).

The equations governing the time dependence of the RyR channel states for the kth Ca²⁺ release subsystem are based on those of Ref. 32, as modified by Ref. 31

(35)

(36)

(37)

(38)

where P_C1,k and P_C2,k are the probabilities of the two RyR closed states for the kth Ca²⁺ release subsystem; m and n are cooperativity parameters; and , , , , , and are rate constants for the transitions between different RyR states (Table 8).

The numerical values of all constants in Eqs. 1–38 are given in Tables 2,3,4,5,6,7,8. Experiments performed at 37°C were adjusted to give appropriate values for room temperature based on the Q₁₀ values for channels, transporters, and pumps given in Ref. 26.

Simulation protocols. Two types of voltage-clamp protocols were simulated to measure activation, inactivation, and recovery from inactivation. Steady-state inactivation was determined with a double-pulse protocol. The first pulse (P1) was a 500-ms depolarization from a holding potential of –75 mV to a voltage between –70 and +50 mV in steps of 5 mV. The membrane potential was then clamped back to –75 mV for 2 ms before the cell was subjected to a second 500-ms pulse (P2) to 0 mV (see Fig. 3A, inset).

View larger version (15K): [in this window] [in a new window]   Fig. 3. Simulated voltage-clamp experiments using a standard double-pulse protocol. A: a 500-ms voltage step (P1) from the holding potential (–75 mV) was applied from –70 to + 50 mV in 10-mV steps, followed by a 500-ms second pulse (P2) to 0 mV (inset). P1 activates the Ca2+ current, which then rapidly inactivates. The magnitude of the peak current elicited by P2 depends on the degree of inactivation at the end of the P1 pulse. B: normalized ICaL (ICaL,norm) versus voltage for both the model (solid line) and experiments [data are from Ref. 55 (), Ref. 52 (), Ref. 53 (), and Ref. 54 ()]. C: steady-state inactivation. The simulated data (solid line) compare well with experimental results [data are from Ref. 53 () and Ref. 13 ()].

Recovery from inactivation was determined using a gapped pulse protocol with a variable interstimulus interval. The first pulse (P1) was a 500-ms depolarization from a holding potential of –75 to 0 mV. This was followed by a second 100-ms depolarization to 0 mV, with the interstimulus interval (when the cell was repolarized to –75 mV) ranging from 2 to 502 ms, in 10-ms intervals.

	RESULTS

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX A: DIFFERENTIAL... APPENDIX B: CALCIUM CURRENTS APPENDIX C: CALCIUM FLUXES REFERENCES

 
L-type Ca²⁺ current. We used the model myocyte to simulate voltage-clamp experiments on isolated ferret ventricular myocytes. The model was subjected to a standard two-pulse protocol (Fig. 3). P1 activated the Ca²⁺ current, and P2 gave an indication of the steady-state inactivation. The simulated results mimic typical voltage-clamp I_CaL experiments (54, 69, 71, 78). The inactivation of I_CaL was biexponential when depolarized to voltages between about –15 and 30 mV. Normalized simulated and experimental current-voltage relationships from ferret myocytes are shown in Fig. 3B. The experimental and modeled results are similar, with a threshold at approximately –30 mV and peak inward current between 0 and +5 mV.

Figure 3C shows simulated (solid line) and experimental values (symbols) of the steady-state inactivation, plotted as a function of the P1 voltage. Steady-state inactivation was calculated by measuring the ratio of the peak current elicited by P2 to the maximum P2 peak current (i.e., peak P2 when P1 was –75 mV). In both the experiment and simulations, inactivation decreased when P1 was positive to about –40 mV and reached a minimum at about 0 mV. Positive to 0 mV inactivation decreased, giving the characteristic U-type inactivation function of ferret I_CaL (53, 54). Although this is a defining characteristic of I_CaL, not all models reproduce this feature.

At potentials positive to about –20 mV, the L-type Ca²⁺ channel inactivates with a biexponential time course. Figure 4A shows simulated and experimental data for the fast and slow inactivation rate constants 1/₁ and 1/₂ (51). The value of the simulated fast rate constant (open squares) is in good agreement with the experimental values (filled squares) in the voltage range of –15 to +40 mV (Fig. 4A). The simulated slow rate constant (open circles) deviates from the experimental values (filled circles) near threshold (Fig. 4A).

View larger version (14K): [in this window] [in a new window]   Fig. 4. A: voltage dependence of ICaL rate constants of inactivation (1/1 and 1/2). ICaL decays in a biexponential manner in response to a depolarizing step. Simulated () and experimental data () for the fast rate constant 1/1 and simulated () and experimental data () for the slow rate constant 1/2 are shown. B: ratio of amplitudes of the fast (A1) and slow (A2) components of ICaL versus voltage from simulated () and experimental data (). The experimental results are from Qu (51).

The voltage dependence of the ratio of amplitudes of the fast and slow components (A₁/A₂) is shown in Fig. 4B. The simulated ratios (open squares) are qualitatively similar to the experimental ones (filled squares), although they show some deviation in amplitude, particularly when one time constant dominates (Fig. 4B).

Recovery from inactivation. We examined the simulated recovery from inactivation with the use of a standard variable interval gapped pulse protocol. An initial pulse to 0 mV was followed by a second pulse to 0 mV after a variable interval, t. The degree of recovery from inactivation was calculated by comparing the peak current elicited by the P2 pulse to that elicited by the P1 pulse. Figure 5A shows the simulated recovery of I_CaL from inactivation, which compares well with experimentally obtained data (53, 54). Figure 5B shows peak simulated I_CaL as a function of interstimulus interval. To fit these data, we used three different functions that have been used to analyze experimental data. These functions were a simple exponential, an exponential function to the 1.5 power and an exponential function squared.

View larger version (21K): [in this window] [in a new window]   Fig. 5. Rate of recovery from inactivation of simulated ICaL. A: simulated ICaL in response to a gapped double-pulse protocol. A 500-ms depolarization to 0 mV was followed by a second 100-ms depolarization to 0 mV after a variable interstimulus interval (t) during which the cell was held at –80 mV. The interstimulus interval changed from 2 to 502 ms, in 25-ms increments. B: simulated peak ICaL from the protocol in A (with 5-ms increments in t) plotted against the interstimulus interval (x). Three different functions were fit to the data: a single exponential f1 = [1 – exp(–t/rec,1)], where rec is the time constant of recovery (thick solid line); an exponential raised to the 1.5 power f2 = [1 – exp(–t/rec,2)]1.5 (thin solid line); and an exponential squared f3 = [1 – exp(–t/rec,3)]2 (dotted line).

All three functions gave a good fit to the modeled recovery curve. Fitting our recovery by a simple exponent gave a recovery time constant of _rec,1 = 98.7 ms. Using the exponential to the 1.5 power gave a recovery time constant of _rec,2 = 77.4 ms, which compares well to the experimental value of _rec,2 = 75 ms (13). The time constant of recovery for the squared exponential was _rec,3 = 66.6 ms, which coincides well with the experimental value of _rec,3 = 68.4 ± 15.1 ms (54).

[Ca²⁺]_i and calcium sparks. Figure 6A shows the modeled [Ca²⁺]_i transients in response to P1 from the same standard voltage-clamp protocol used in Fig. 3. The peak value of the simulated [Ca²⁺]_i elicited by the first pulse plotted as a function of voltage is shown in Fig. 6B. The bell-shaped relationship shown by the model is similar to the voltage dependence of peak [Ca²⁺]_i of experimentally observed cardiac myocytes from various species (1, 9, 15, 19, 45, 74). The bell-shaped relationship between [Ca²⁺]_i and voltage reflects graded release.

View larger version (15K): [in this window] [in a new window]   Fig. 6. Simulated Ca2+ transients in response to the double-pulse protocol shown in Fig. 3. A: [Ca2+]i. B: peak value of [Ca2+]i ([Ca2+]i,max) in P1 as function of P1 voltage.

The calcium concentrations [Ca²⁺]_ss,k (k = 1–8) in the eight subcellular volumes have independent time courses. The number of subsystems in which CICR is initiated varies with the subsystem type and the transmembrane voltage. Figure 7 shows the time course of [Ca²⁺]_ss,k for all the Ca²⁺ release subsystems as well as the general [Ca²⁺]_i for various P1 voltages. When P1 was negative to about –25 mV (i.e., below threshold for I_CaL activation), there was almost no change in either [Ca²⁺]_ss,k or [Ca²⁺]_i. If P1 was more depolarized, e.g., –20 mV, CICR occurred in one subsystem, leading to a single Ca²⁺ spark in one subspace volume. This rapid local increase of [Ca²⁺]_ss,k was followed by a slow increase in [Ca²⁺]_i to 0.14 µM. When P1 was –15 mV, there were four Ca²⁺ sparks and a fairly slow increase in [Ca²⁺]_i to 0.33 µM (Fig. 7A). All eight calcium release subsystems produce Ca²⁺ sparks, i.e., the maximum number of sparks possible, when P1 was between 0 and +15 mV (Fig. 7B). This was accompanied by the fastest increase in [Ca²⁺]_i, which reached a peak of 0.76 µM within 25 ms. When P1 was further depolarized to a value such as +20 mV, there was a decrease in the number of Ca²⁺ subsystems that release Ca²⁺ from the SR, to seven. There was a concurrent decrease in rate and the maximum value of [Ca²⁺]_i (0.72 µM at 28 ms). When P1 was +35 mV, six of the eight potential populations of release units produced sparks and [Ca²⁺]_i peaked at 0.59 µM (Fig. 7C). As the voltage of P1 was further depolarized, fewer sparks were initiated, until at +50 mV only two sparks were seen, and there was only a small slow increase in [Ca²⁺]_i (Fig. 7D). Our model therefore shows gradation of Ca²⁺ release.

View larger version (28K): [in this window] [in a new window]   Fig. 7. Relationship between subspace CICR release and bulk [Ca2+]i transients during depolarization. The number of subsystems that exhibit Ca2+ transients (i.e., sparks) in response to a depolarization changes with voltage. The subspace Ca2+ concentration ([Ca2+]ss) for all eight subsystems is shown. A: when P1 is –15 mV, there are few sparks and [Ca2+]i is delayed. B: peak [Ca2+]i occurs rapidly near 0 mV. C: when P1 is +35 mV, [Ca2+]i is delayed and there are fewer sparks. D: when P1 is +50 mV, there is a slow rise in [Ca2+]i and calcium transients occur in only a few subsystems.

An interesting result of our model is the relationship between the P1 voltage and the time over which sparks occur. At voltages just positive to the threshold of activation of I_CaL (P1 = –20 to –10 mV), Ca²⁺ sparks are initiated over relatively wide time window, up to 40 ms after the beginning of voltage pulse (Fig. 7A). When P1 is near the maximum for I_CaL, most sparks are initiated earlier. As P1 is depolarized further (35–50 mV), sparks occur in a less concentrated fashion (Fig. 7C). A similar relationship between depolarization voltage and the time to initiation of Ca²⁺ sparks was observed experimentally in guinea pig ventricular myocytes (38). In this experiment, Ca²⁺ sparks were observed over a wide time window both at relatively small voltage steps near the threshold of Ca²⁺ release from the SR and at relatively large voltage steps (i.e., P1 = 30 mV). However, when P1 is in the range 0–10 mV, most Ca²⁺ sparks are concentrated near the beginning of the depolarization.

The number of Ca²⁺ release subsystems that produce Ca²⁺ sparks varies with voltage. The percentage of available Ca²⁺ release systems in which sparks occur (the number of release systems that spark multiplied by the number of those systems present) produces the net Ca²⁺ transient. Despite the relatively small number of Ca²⁺ release subsystems used in this model, the relationship between number of sparks and voltage is bell shaped, which is similar to the experimentally obtained bell-shape dependence of local [Ca²⁺]_i transients on voltage (38). These data suggest that the inhomogeneity of Ca²⁺ release subsystems is responsible for gradedness of Ca²⁺ release: Ca²⁺ release subsystems with different properties distributed nonuniformly in a cardiac cell produce nonuniformity in Ca²⁺ release and [Ca²⁺]_i. Spatial nonuniformities between [Ca²⁺]_i and the number of Ca²⁺ sparks have been observed in cardiac cells (16, 74). The role of Ca²⁺ release subsystem heterogeneity in generating graded Ca²⁺ release is clearly shown by the lack of graded release when all subsystems are homogeneous with 40 RyRs and 4 DHPRs in each subsystem (Fig. 8, squares). When only four types of release subsystems (with 2, 4, 6, and 8 DHPR and 40 RyRs in each subsystem; Fig. 8, triangles) are modeled, the relationship between voltage and [Ca²⁺]_i is less bell shaped than with eight subsystems.

View larger version (14K): [in this window] [in a new window]   Fig. 8. Effect of changing the number of subsystems on graded release. Peak normalized [Ca2+]i is plotted as a function of voltage for three different model configurations. The homogeneous case, with 4 L-type Ca2+ channels, showed no gradation of release (). When there were 4 release subsystem types, with either 2, 4, 6, or 8 L-type Ca2+ channels, mild gradation was found (). With 8 subsystems, gradation of Ca2+ release was smoother ().

Factors affecting graded release. The sensitivity of graded release to the total number of subsystems modeled in this study suggested that graded release may also be sensitive to other factors. We examined the sensitivity of the graded release response to various weighting distributions for the subunit sizes with Ca²⁺ release subsystems containing one to eight DHPRs. These weighting factors followed an exponential (exp_k), uniform (unif_k), or Maxwellian (maxw_k) distribution (weighting factors given in Table 2). The results are shown in Fig. 9. Graded release was smoother and more bell shaped for the exponential and uniform distribution of weighting factors, particularly for negative membrane potentials. This demonstrates that the shape of graded release produced by this model is sensitive to the type of weighting distribution of Ca²⁺ release subsystem size.

View larger version (24K): [in this window] [in a new window]   Fig. 9. Effect of weighting distributions on graded release. Normalized ICaL amplitudes () and peak [Ca2+]i () are shown as a function of voltage. A–D: calculations for binomial, uniform, exponential, and Maxwellian distributions of weights for Ca2+ release subsystems, respectively.

Other factors have been suggested to play an important role in graded release. One of these factors is RyR adaptation. The role of adaptation was evaluated by setting the rate constant to zero, which corresponds to removal of adaptation (Fig. 10). In agreement with the study of Jafri et al. (31), this manipulation had little effect on the maximal Ca²⁺ release transient and time course (Fig. 10, B–D). However, adaptation did alter the shape of the graded release response. This suggests that although adaptation does not give rise to graded release in our model, it may still play a significant role in modulating the properties of graded release.

View larger version (12K): [in this window] [in a new window]   Fig. 10. Effect of RyR adaptation on graded release. Simulations were run with either default RyRs (with adaptation) or RyRs with no adaptation. A 500-ms depolarizing step from –75 to 0 mV was applied to each model. A: normalized ICaL amplitudes (squares) and peak [Ca2+]i (triangles) as a function of voltage. Filled symbols show control adaptation data [RyR PO1,k – PC2,k rate constant and open symbols are data calculated in the absence of adaptation . There was no change in the voltage dependence of peak ICaL, and squares superimpose. B: time dependence of [Ca2+]i in response to step depolarization. Thick solid lines, control data; thin solid lines, data obtained with no RyR adaptation. C: time dependence of Ca2+ concentration in the JSR facing a sarcolemmal region with 4 L-type Ca2+ channels present ([Ca2+]JSR,4) in response to step depolarization. Thick solid lines, control data; thin solid lines, data from RyRs without adaptation. D: time dependence of the open probability of RyRs in the 4 L-type Ca2+ channel subsystem (PRyRopen,4) in response to step depolarization. Thick solid lines, control data; thin solid lines, nonadapting RyR data.

Effect of SR load on graded release. The amount of Ca²⁺ in the SR is an important and sensitive determinant of Ca²⁺ release (6). We increased and decreased the Ca²⁺ load in this model by a factor of 2. Figure 11A shows that the total amount of Ca²⁺ released was approximately proportional to SR load. The loaded SR increased [Ca²⁺]_ss in the subspace volume and prolonged the Ca²⁺ spark, which in turn kept the RyR in an open state for a longer period of time. Interestingly, release became less graded as the SR load increased, with maximal Ca²⁺ release occurring with a much smaller Ca²⁺ entry. Thus it seems that increased SR load makes the RyR more likely to open. This occurs without making the opening rate constants for the RyR explicitly dependent on SR Ca²⁺ load, unlike the model of Sobie et al. (64). In our model, the increased sensitivity comes from the regenerative nature of Ca²⁺ release in the individual subsystem. Within each subsystem when Ca²⁺ in the junctional space reaches a critical level, Ca²⁺ entry through RyRs causes regenerative activation of further RyRs. When the SR is loaded, the driving force for Ca²⁺ diffusion is greatly increased. As a consequence, a lower value of channel open probability is required to produce the same influx of Ca²⁺ into the subspace. Therefore, there is a lower threshold for triggering all-or-none Ca²⁺ release in a subsystem when the SR is Ca²⁺ loaded.

View larger version (19K): [in this window] [in a new window]   Fig. 11. Modulation of SR Ca2+ release by SR load and RyR density. A: effect of SR load on CICR. Calculations were made with 50%, 100%, and 200% of control SR load. Filled symbols show the relationship between ICaL and voltage, normalized to the maximum value with 100% of SR load (control). There is little difference in this relationship with changes in SR loading. Open symbols show the relationship between [Ca2+]i and voltage, normalized to the peak [Ca2+]i obtained with 100% SR loading (control, squares). Triangles, data from 200% SR loading; circles, data from 50% SR loading. B: effect of varying the number of RyRs per release subsystem. RyR conductance was calculated as 50%, 100%, and 200% of control. The normalized ICaL amplitudes (filled symbols) are little effected by changes in RyR density. Open symbols show the relationship between [Ca2+]i and voltage for control (100% RyR density, squares) and high (200% density, triangles) and low (50% density, circles) RyR density.

Altering the total number of RyRs had similar effects on changing of the SR Ca²⁺ load. Increasing the number of RyRs increased Ca²⁺ release (Fig. 11B, open triangles), increased the rate of inactivation of I_CaL, and resulted in graded release that was more abrupt and had near-maximal release with a much smaller trigger Ca²⁺. However, changing the relative weighting so that there were more RyRs in larger release units and fewer RyRs in smaller units caused little change in the graded response (simulations not shown).

Reducing Ca²⁺ influx via the L-type Ca²⁺ channel by 50% resulted in a graded response with a diminished and delayed transient (Fig. 12). Consistent with a decrease in Ca²⁺ release, the inactivation rate of the channel was also slowed. The slowing of the inactivation rate and the decreased intracellular Ca²⁺ release are consistent with changes seen after block of I_CaL by dihydropyridine antagonists or the reduced I_CaL density accompanying hypertrophic heart failure or myocardial stunning.

View larger version (17K): [in this window] [in a new window]   Fig. 12. Effect of ICaL block on CICR. A: simulated time dependence of ICaL using a standard double-pulse protocol. A 500-ms voltage step (P1) from the holding potential (–75 mV) was applied from –70 to + 50 mV in 10-mV steps, followed by a 500-ms second pulse (P2) to 0 mV. The maximum value of ICaL was reduced by 50% to simulate the action of an ICaL antagonist. B: normalized ICaL as a function of voltage in control () and after ICaL antagonist (). Currents were normalized to the maximum value of control ICaL. Peak [Ca2+]i is also shown as a function of voltage (, control data; , antagonist data). All data were normalized to the maximum value of control [Ca2+]i. C: total Ca2+ flux from the JSR into the subspace for all subsystems under control and antagonist conditions in response to a step depolarization from –75 to 0 mV.

Calcium release during deactivation. Ca²⁺ release from the SR terminates about 40 ms after activation of the L-type Ca²⁺ channel. This may be due to many factors, including stochastic attrition of Ca²⁺ release subunits, Ca²⁺ channel inactivation, and RyR inactivation or adaptation. Sham et al. (61) used the Ca²⁺ channel agonist FPL-64176 to remove Ca²⁺-dependent inactivation, increase the peak current, and slow deactivation of I_CaL. Sham et al. (61) presented data that "... argue against the stochastic attrition (66) as the primary mechanism for terminating Ca²⁺ release, because it predicts an increase in open duration of L-type Ca²⁺ channels would prolong Ca²⁺ release and multiple Ca²⁺ channel reopenings would simply give rise to multiple release events (66)."

Cells treated with FPL-64176 were depolarized to –30, 0, +30, or +60 mV for 50 ms and then hyperpolarized to –120 mV. When the initial depolarization was to –30 or 0 mV, near-maximal Ca²⁺ release was seen. The subsequent hyperpolarization either failed to trigger or triggered only a minimal secondary Ca²⁺ release. However, when the initial depolarization was to +30 or +60 mV (and the primary Ca²⁺ release was submaximal), a substantial "tail" of secondary Ca²⁺ release was seen on hyperpolarization.

We simulated these protocols for –20, 0, and +30 mV (Fig. 13). The simulations reproduced the minimal secondary Ca²⁺ release at 0 mV and the large secondary Ca²⁺ release at +30 mV. This is an important validation of our model. The secondary Ca²⁺ release during the hyperpolarizing step must be from a different population of release units than the primary Ca²⁺ release.

View larger version (15K): [in this window] [in a new window]   Fig. 13. Initiation of a secondary Ca2+ release with inactivation-modified L-type Ca2+ channels by hyperpolarization. These simulations mimic the exposure of rat myocytes to FPL-64176 (FPL) (61). FPL removes fast Ca2+-dependent inactivation and slows deactivation of L-type Ca2+ channels {(V) replaced by '(V) = (V)/[1 + 0.5e–(V – 3.0)/13.0], where V is membrane potential}. A: summed [Ca2+]ss from all subunits for a depolarization to –20 mV followed by hyperpolarization to –120 mV (see inset for protocol). B: corresponding ICaL for A. C: summed [Ca2+]ss from all subunits for a depolarization to 0 mV followed by hyperpolarization to –120 mV. D: corresponding ICaL for C. E: summed [Ca2+]ss from all subunits for a depolarization to +30 mV followed by hyperpolarization to –120 mV. F: corresponding ICaL for E. Note that a separate population of subunits fires when higher single channel currents are elicited on hyperpolarization after a submaximal Ca2+ release during depolarization.

	DISCUSSION

TOP ABSTRACT METHODS RESULTS DISCUSSION APPENDIX A: DIFFERENTIAL... APPENDIX B: CALCIUM CURRENTS APPENDIX C: CALCIUM FLUXES REFERENCES

 
We developed a model of calcium dynamics in the ferret ventricular myocyte that exhibits graded CICR from the SR. It is based on a system of eight distinct Ca²⁺ release subunits, each with slightly different characteristics, and a new Markov model of the L-type Ca²⁺ channel. Our model successfully simulates voltage-clamp experiments on I_CaL and subsequent CICR in ferret ventricular cardiac myocytes both quantitatively and qualitatively.

This model of CICR produces graded release. The mechanism by which it does so is based on diad inhomogeneity. The total number of Ca²⁺ channels per diad is not uniform, but there is a binomial distribution, with four Ca²⁺ channels being the average and most common. For a small Ca²⁺ influx, the subspace Ca²⁺ will rise only slightly in diads with only one or two Ca²⁺ channels but will rise much higher in diads with seven or eight Ca²⁺ channels. Thus for small I_CaL relatively few subunits reach the Ca²⁺ release threshold, and little Ca²⁺ release is expected. When Ca²⁺ influx is higher, diads with fewer channels will reach the threshold for Ca²⁺ release. With large I_CaL virtually all available units release. This model is based on biological variability.

The magnitude of the [Ca²⁺]_i transient depends on the magnitude of I_CaL (for reviews, see Refs. 11 and 73). Our model exhibits a nonlinear relationship between [Ca²⁺]_i transients and I_CaL, which has also been seen experimentally (12, 14). Figure 14 shows that there are bell-shaped relationships among normalized I_CaL, normalized [Ca²⁺]_i, normalized number of Ca²⁺ sparks (n_sparks; or normalized number of Ca²⁺ release subsystems that produce Ca²⁺ sparks), and voltage. The curves for [Ca²⁺]_i and n_sparks are similar, whereas the relationship for I_CaL is comparatively narrow. This indicates that [Ca²⁺] in the general cytosol is determined mainly by the magnitude of CICR from the SR and not I_CaL. Our model shows a degree of relative independence among I_CaL, Ca²⁺ dynamics, and Ca²⁺ release from the SR, an observation that has also been determined experimentally (61).

View larger version (16K): [in this window] [in a new window]   Fig. 14. Normalized ICaL amplitudes (), peak [Ca2+]i (), and number of Ca2+ release subsystems that produce Ca2+ sparks (nsparks, ) as a function of voltage.

The eight different Ca²⁺ subsystems in our model results in an inhomogeneity in local Ca²⁺ release, which results in graded CICR, and the staggered appearance of Ca²⁺ sparks. This is an experimentally observed feature (38) and was not explained by any previous model. Our model also explains the relationship between membrane voltage and the number of Ca²⁺ sparks initiated (38).

Two previously published models have simulated graded CICR from the SR: the Luo-Rudy model (39) and the Rice et al. model (57). The Luo-Rudy model does not consider the local nature of Ca²⁺ release from the SR and therefore has no information on Ca²⁺ sparks. The Rice et al. model (57) includes some features of