## Abstract

We have developed a model of Ca^{2+} handling in ferret ventricular myocytes. This model includes a novel L-type Ca^{2+} channel, detailed intracellular Ca^{2+} movements, and graded Ca^{2+}-induced Ca^{2+} release (CICR). The model successfully reproduces data from voltage-clamp experiments, including voltage- and time-dependent changes in intracellular Ca^{2+} concentration ([Ca^{2+}]_{i}), L-type Ca^{2+} channel current (*I*_{CaL}) inactivation and recovery kinetics, and Ca^{2+} 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^{2+} 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^{2+} channels distributed binomially. This model exhibits graded CICR and provides a quantitative description of Ca^{2+} dynamics not requiring Monte-Carlo simulations. Activation of RyRs and release of Ca^{2+} from the SR depend critically on Ca^{2+} entry through L-type Ca^{2+} channels. In turn, Ca^{2+} channel inactivation is critically dependent on the release of stored intracellular Ca^{2+}. Inactivation of *I*_{CaL} depends on both transmembrane voltage and local [Ca^{2+}]_{i} near the channel, which results in distinctive inactivation properties. The molecular mechanisms underlying many *I*_{CaL} gating properties are unclear, but [Ca^{2+}]_{i} dynamics clearly play a fundamental role.

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

calcium is of vital importance in cells as both a transmembrane electrical signaling ion and a second messenger (3–5, 7, 36, 50, 59). The diverse range of Ca^{2+}-mediated effects is particularly evident in cardiac myocytes, where the initial transmembrane flux of ions through L-type Ca^{2+} channels contributes to the action potential (6, 25, 40, 46). The entering Ca^{2+} then initiates Ca^{2+}-induced Ca^{2+} release (CICR) from the sarcoplasmic reticulum (SR), which gives rise to the intracellular Ca^{2+} concentration ([Ca^{2+}]_{i}) transient, which in turn affects muscle contraction, transmitter release, gene transcription, channel modulation, etc. (5, 10, 18, 36, 59, 73). [Ca^{2+}]_{i} cycling is a fundamental part of normal cellular function, so any malfunction or perturbation of the systems that handle Ca^{2+} can result in heart failure or fibrillation (21, 33, 44, 45, 49, 58, 60, 67, 69, 75).

Recently, CICR has been observed at the subcellular level in the form of Ca^{2+} “sparks,” which are extremely large, extremely localized increases in [Ca^{2+}] (18, 19, 29, 34, 35). Spontaneously occurring sparks have been observed. They can trigger a subcellular Ca^{2+} wave that can propagate throughout the cell (17, 30, 35). Ca^{2+} 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^{2+} 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^{2+}] cycling is the graded release of Ca^{2+} 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^{2+} 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^{2+} handling that incorporates a novel Markov model for L-type Ca^{2+} channels and graded SR Ca^{2+} 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^{2+} channels, and the junctional SR. The SR surface exposed to the subsystem has a cluster of RyRs that release Ca^{2+} into the subsystem volume. The subsystems are distributed throughout the model cell, with the number of L-type Ca^{2+} 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^{2+} entering the subspace via the L-type Ca^{2+} channels is sufficiently large. We also developed a new Markov model for the L-type Ca^{2+} channel. This L-type Ca^{2+} 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^{2+}- 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^{2+} sparks observed (as determined by confocal scanning microscopy experiments).

## METHODS

*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^{2+} channels, which allow Ca^{2+} to enter the subspace volume (V_{ss}) and increase the subspace Ca^{2+} concentration ([Ca^{2+}]_{ss}). The entering calcium (*I*_{CaL}) triggers CICR from the junctional SR (*J*_{rel}), which induces further CICR and inactivates *I*_{CaL} (65). Ca^{2+} in the subspace diffuses from the subsystem to the general cytosol (*J*_{xfer}), where it contributes to [Ca^{2+}]_{i}. Calcium leaves the general cytosol by being pumped back into the SR by SR Ca^{2+}-ATPase (*J*_{up}), pumped across the cell membrane by the sarcolemmal Ca^{2+} pump [Ca^{2+} pump current (*I*_{pCa})], or extruded from the cell by the Na^{+}/Ca^{2+} exchange mechanism [Na^{+}/Ca^{2+} exchange current (*I*_{NaCa})]. Calcium buffers present in the general cytosol also remove Ca^{2+} from free cytosolic activity: Ca^{2+} 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^{2+} due to the presence of calsequestrin. Calcium is assumed to be in instantaneous equilibrium with calsequestrin. A small leak current takes Ca^{2+} directly from the network SR to the general cytosol (*J*_{leak}), and the background Ca^{2+} current (*I*_{Cab}) transfers Ca^{2+} from the extracellular space to the main cytosol.

*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^{2+} flux through the L-type Ca^{2+} channel enters, CICR is initiated, and SR Ca^{2+} is released. Each subsystem has between one and eight L-type Ca^{2+} channels and a cluster of RyRs. Figure 1*B* shows three separate Ca^{2+} release subsystems, Ca_{ss5}, Ca_{ss1}, and Ca_{ss8}, with five, one, and eight L-type Ca^{2+} 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^{2+} channels per RyR cluster. The weighting functions for L-type Ca^{2+} channel distributions (bin* _{k}*, unif

*, exp*

_{k}*, and maxw*

_{k}*) are described in Table 2.*

_{k}

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 × 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^{2+} 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 × 10^{–9} μl (360 nm × 360 nm × 12 nm × 1,250). There is a linear relationship between V_{ss,}* _{k}* and the number of L-type Ca

^{2+}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

^{2+}channels).

For the subsystems with *k* L-type Ca^{2+} channels, the L-type Ca^{2+} current is *I*_{CaL,k}, the Ca^{2+} concentration in V_{ss,k} is [Ca^{2+}]_{ss,}* _{k}*, the CICR flux is

*J*

_{rel,}

*, and the diffusion of Ca*

_{k}^{2+}from the subspace volume to the main cytosol is

*J*

_{xfer,k}. V

_{ss,}

*depends linearly on the number of L-type Ca*

_{k}^{2+}channels present. The time constant of the transfer of Ca

^{2+}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

^{2+}efflux from the SR,

*J*

_{rel,}

*, is different for each subsystem.*

_{k}

The total Ca^{2+} flux depends on the sum of the activity in all subsystems with *k* L-type Ca^{2+} 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

^{2+}] between the junctional SR ([Ca

^{2+}]

_{JSR,}

*) and the subspace volume ([Ca*

_{k}^{2+}]

_{ss,}

*), the sum of probabilities for opening of the RyRs (*

_{k}*P*

_{O1,}

*and*

_{k}*P*

_{O2,}

*), and the maximum permeability of RyR channels (*

_{k}*v*

_{1}) (3)

The magnitude of the Ca^{2+} flux from the network to junctional SR (*J*_{tr,}* _{k}*) depends on the differences in [Ca

^{2+}] (i.e., [Ca

^{2+}]

_{NSR}and [Ca

^{2+}]

_{JSR,}

*) and the time constant of diffusion (τ*

_{k}_{tr}) (4)

Ca^{2+} flux from the subspace volume to the cytoplasm (*J*_{xfer,k}) depends on the relative concentrations of calcium in the subspace ([Ca^{2+}]_{ss,k}) and the cytosol ([Ca^{2+}]_{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^{2+} flux to the SR (*J*_{up}) depends on [Ca^{2+}]_{i}, the maximum SR Ca^{2+}-ATPase pump rate (*v*_{3}), and the half-saturation constant for SR Ca^{2+}-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^{2+} bound, [HTRPNCa] is the concentration of high-affinity troponin-binding sites with Ca^{2+} bound, is the Ca^{2+} on rate for troponin high-affinity sites, is the Ca^{2+} off rate for troponin high-affinity sites, is the Ca^{2+} on rate for low-affinity troponin sites, and is the Ca^{2+} off rate for low-affinity troponin sites (Table 4). Calmodulin, the other cytosolic Ca^{2+} buffer, is assumed to be in instantaneous equilibrium with [Ca^{2+}] (see *Eqs. 15, 16, 19*, and *20*).

The leak of Ca^{2+} from the SR to the cytosol (*J*_{leak}) depends on the relative concentrations of Ca^{2+} and the leak rate constant (*v*_{2}) (8)

*L-type Ca ^{2}*

^{+}

*channels.*Modeling the inactivation behavior of the L-type Ca

^{2+}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

^{2+}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

^{2+}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

^{2+}-dependent inactivation and slow Ca

^{2+}-insensitive inactivation. Fast calcium- and calmodulin-sensitive inactivation may be similar to N-type inactivation (48) and slow Ca

^{2+}-independent inactivation may be similar to C-type inactivation (77). This formulation can reproduce the behavior of the Ca

^{2+}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^{2+} channel with four closed states (C_{1}, C_{2}, C_{3}, and C_{4}), one open state (O), and three inactivated states (I_{1}, I_{2}, and I_{3}). 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^{2+}-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_{2} transition represents the channel undergoing “voltage-dependent” inactivation. Fast Ca^{2+}-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^{2+} to bind to “activate” it (20). The O → I_{1} transition corresponds to Ca^{2+}-dependent inactivation, which is controlled by the Ca^{2+}-dependent rate constant γ* _{k}*. I

_{3}represents a channel that is both voltage and Ca

^{2+}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

_{4}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.

Each Ca^{2+} subsystem contains a different number of L-type Ca^{2+} channels and different subspace Ca^{2+} concentrations, [Ca^{2+}]_{ss,}* _{k}*. The calcium-induced inactivation of the channel, and therefore the behavior of the L-type Ca

^{2+}channel, will be different in each subsystem. The total

*I*

_{CaL}is the sum of all the subsystem Ca

^{2+}currents (9) where

*N*is the number of calcium release subsystems and

*I*

_{CaL,}

*is*

_{k}*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

*is the probability of the channel being open.*

_{k}*Other Ca ^{2}*

^{+}

*currents.*The model includes three other Ca

^{2+}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

^{2+}half-saturation constant.

*I*

_{Cab}is modeled as a simple linear ohmic conductor (12) (13) where

*G*

_{Cab}is background Ca

^{2+}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^{2+}, 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^{2+} half-saturation constant, *K*_{sat} is the Na^{+}/Ca^{2+} 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^{2+}] 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

^{2+}half-saturation constant for calmodulin, [CSQN]

_{tot}is the network SR calsequestrin concentration, and is the Ca

^{2+}half-saturation constant for calsequestrin. Note that the fluxes and variables with a

*k*index in

*Eqs. 15–17*refer to the

*k*th Ca

^{2+}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^{2+} binding to calsequestrin and calmodulin is considered to be instantaneous (31). In contrast, Ca^{2+} 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^{2+}]_{ss,}* _{k}*, so the state of the L-type Ca

^{2+}channels will vary between subsystems. The time course of the channel for the conducting state (O

*), the closed states (C*

_{k}_{1,}

*–C*

_{k}_{4,}

*), and the inactivated states (I*

_{k}_{1,}

*–I*

_{k}_{3,}

*) for the*

_{k}*k*th Ca

^{2+}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

^{2+}-induced inactivation and

*K*

_{pc,half}is the half-saturation constant for Ca

^{2+}-induced inactivation (listed in Table 7).

The equations governing the time dependence of the RyR channel states for the *k*th Ca^{2+} 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,}

*are the probabilities of the two RyR closed states for the*

_{k}*k*th Ca

^{2+}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_{10} 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. 3*A*, *inset*).

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

*L-type Ca ^{2}*

^{+}

*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

^{2+}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. 3

*B*. The experimental and modeled results are similar, with a threshold at approximately –30 mV and peak inward current between 0 and +5 mV.

Figure 3*C* 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^{2+} channel inactivates with a biexponential time course. Figure 4*A* shows simulated and experimental data for the fast and slow inactivation rate constants 1/τ_{1} and 1/τ_{2} (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. 4*A*). The simulated slow rate constant (open circles) deviates from the experimental values (filled circles) near threshold (Fig. 4*A*).

The voltage dependence of the ratio of amplitudes of the fast and slow components (A_{1}/A_{2}) is shown in Fig. 4*B*. 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. 4*B*).

*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 5*A* shows the simulated recovery of *I*_{CaL} from inactivation, which compares well with experimentally obtained data (53, 54). Figure 5*B* 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.

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 ^{2}*

^{+}

*]*Figure 6

_{i}and calcium sparks.*A*shows the modeled [Ca

^{2+}]

_{i}transients in response to P1 from the same standard voltage-clamp protocol used in Fig. 3. The peak value of the simulated [Ca

^{2+}]

_{i}elicited by the first pulse plotted as a function of voltage is shown in Fig. 6

*B*. The bell-shaped relationship shown by the model is similar to the voltage dependence of peak [Ca

^{2+}]

_{i}of experimentally observed cardiac myocytes from various species (1, 9, 15, 19, 45, 74). The bell-shaped relationship between [Ca

^{2+}]

_{i}and voltage reflects graded release.

The calcium concentrations [Ca^{2+}]_{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

^{2+}]

_{ss,k}for all the Ca

^{2+}release subsystems as well as the general [Ca

^{2+}]

_{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

^{2+}]

_{ss,k}or [Ca

^{2+}]

_{i}. If P1 was more depolarized, e.g., –20 mV, CICR occurred in one subsystem, leading to a single Ca

^{2+}spark in one subspace volume. This rapid local increase of [Ca

^{2+}]

_{ss,}

*was followed by a slow increase in [Ca*

_{k}^{2+}]

_{i}to 0.14 μM. When P1 was –15 mV, there were four Ca

^{2+}sparks and a fairly slow increase in [Ca

^{2+}]

_{i}to 0.33 μM (Fig. 7

*A*). All eight calcium release subsystems produce Ca

^{2+}sparks, i.e., the maximum number of sparks possible, when P1 was between 0 and +15 mV (Fig. 7

*B*). This was accompanied by the fastest increase in [Ca

^{2+}]

_{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

^{2+}subsystems that release Ca

^{2+}from the SR, to seven. There was a concurrent decrease in rate and the maximum value of [Ca

^{2+}]

_{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

^{2+}]

_{i}peaked at 0.59 μM (Fig. 7

*C*). 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

^{2+}]

_{i}(Fig. 7

*D*). Our model therefore shows gradation of Ca

^{2+}release.

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^{2+} sparks are initiated over relatively wide time window, up to 40 ms after the beginning of voltage pulse (Fig. 7*A*). 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. 7*C*). A similar relationship between depolarization voltage and the time to initiation of Ca^{2+} sparks was observed experimentally in guinea pig ventricular myocytes (38). In this experiment, Ca^{2+} sparks were observed over a wide time window both at relatively small voltage steps near the threshold of Ca^{2+} 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^{2+} sparks are concentrated near the beginning of the depolarization.

The number of Ca^{2+} release subsystems that produce Ca^{2+} sparks varies with voltage. The percentage of available Ca^{2+} 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^{2+} transient. Despite the relatively small number of Ca^{2+} 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^{2+}]_{i} transients on voltage (38). These data suggest that the inhomogeneity of Ca^{2+} release subsystems is responsible for gradedness of Ca^{2+} release: Ca^{2+} release subsystems with different properties distributed nonuniformly in a cardiac cell produce nonuniformity in Ca^{2+} release and [Ca^{2+}]_{i}. Spatial nonuniformities between [Ca^{2+}]_{i} and the number of Ca^{2+} sparks have been observed in cardiac cells (16, 74). The role of Ca^{2+} release subsystem heterogeneity in generating graded Ca^{2+} 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^{2+}]_{i} is less bell shaped than with eight subsystems.

*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^{2+} release subsystems containing one to eight DHPRs. These weighting factors followed an exponential (exp* _{k}*), uniform (unif

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

_{k}^{2+}release subsystem size.

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^{2+} 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.

*Effect of SR load on graded release.* The amount of Ca^{2+} in the SR is an important and sensitive determinant of Ca^{2+} release (6). We increased and decreased the Ca^{2+} load in this model by a factor of 2. Figure 11*A* shows that the total amount of Ca^{2+} released was approximately proportional to SR load. The loaded SR increased [Ca^{2+}]_{ss} in the subspace volume and prolonged the Ca^{2+} 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^{2+} release occurring with a much smaller Ca^{2+} 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^{2+} load, unlike the model of Sobie et al. (64). In our model, the increased sensitivity comes from the regenerative nature of Ca^{2+} release in the individual subsystem. Within each subsystem when Ca^{2+} in the junctional space reaches a critical level, Ca^{2+} entry through RyRs causes regenerative activation of further RyRs. When the SR is loaded, the driving force for Ca^{2+} diffusion is greatly increased. As a consequence, a lower value of channel open probability is required to produce the same influx of Ca^{2+} into the subspace. Therefore, there is a lower threshold for triggering all-or-none Ca^{2+} release in a subsystem when the SR is Ca^{2+} loaded.

Altering the total number of RyRs had similar effects on changing of the SR Ca^{2+} load. Increasing the number of RyRs increased Ca^{2+} release (Fig. 11*B*, 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^{2+}. 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^{2+} influx via the L-type Ca^{2+} channel by 50% resulted in a graded response with a diminished and delayed transient (Fig. 12). Consistent with a decrease in Ca^{2+} release, the inactivation rate of the channel was also slowed. The slowing of the inactivation rate and the decreased intracellular Ca^{2+} 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.

*Calcium release during deactivation.* Ca^{2+} release from the SR terminates about 40 ms after activation of the L-type Ca^{2+} channel. This may be due to many factors, including stochastic attrition of Ca^{2+} release subunits, Ca^{2+} channel inactivation, and RyR inactivation or adaptation. Sham et al. (61) used the Ca^{2+} channel agonist FPL-64176 to remove Ca^{2+}-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^{2+} release, because it predicts an increase in open duration of L-type Ca^{2+} channels would prolong Ca^{2+} release and multiple Ca^{2+} 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^{2+} release was seen. The subsequent hyperpolarization either failed to trigger or triggered only a minimal secondary Ca^{2+} release. However, when the initial depolarization was to +30 or +60 mV (and the primary Ca^{2+} release was submaximal), a substantial “tail” of secondary Ca^{2+} release was seen on hyperpolarization.

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

## DISCUSSION

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^{2+} release subunits, each with slightly different characteristics, and a new Markov model of the L-type Ca^{2+} 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^{2+} channels per diad is not uniform, but there is a binomial distribution, with four Ca^{2+} channels being the average and most common. For a small Ca^{2+} influx, the subspace Ca^{2+} will rise only slightly in diads with only one or two Ca^{2+} channels but will rise much higher in diads with seven or eight Ca^{2+} channels. Thus for small *I*_{CaL} relatively few subunits reach the Ca^{2+} release threshold, and little Ca^{2+} release is expected. When Ca^{2+} influx is higher, diads with fewer channels will reach the threshold for Ca^{2+} release. With large *I*_{CaL} virtually all available units release. This model is based on biological variability.

The magnitude of the [Ca^{2+}]_{i} transient depends on the magnitude of *I*_{CaL} (for reviews, see Refs. 11 and 73). Our model exhibits a nonlinear relationship between [Ca^{2+}]_{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^{2+}]_{i}, normalized number of Ca^{2+} sparks (*n*_{sparks}; or normalized number of Ca^{2+} release subsystems that produce Ca^{2+} sparks), and voltage. The curves for [Ca^{2+}]_{i} and *n*_{sparks} are similar, whereas the relationship for *I*_{CaL} is comparatively narrow. This indicates that [Ca^{2+}] 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^{2+} dynamics, and Ca^{2+} release from the SR, an observation that has also been determined experimentally (61).

The eight different Ca^{2+} subsystems in our model results in an inhomogeneity in local Ca^{2+} release, which results in graded CICR, and the staggered appearance of Ca^{2+} 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^{2+} 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^{2+} release from the SR and therefore has no information on Ca^{2+} sparks. The Rice et al. model (57) includes some features of local Ca^{2+} release events and simulates Ca^{2+} sparks after rapid Ca^{2+} efflux from the SR. The Rice et al. model (57) has a stochastically opening L-type Ca^{2+} channel that produces graded SR Ca^{2+} release. In their model, an increase in the open probability of the L-type Ca^{2+} channel is accompanied by an increase in CICR from the SR. However, the experimental results of Sham et al. (61) showed that use of an agonist to increase the open probability of the L-type Ca^{2+} channel neither prolongs CICR nor triggers secondary Ca^{2+} sparks.

What sort of advantages would spatial heterogeneity give a cell relative to a purely stochastic graded release mechanism? We can speculate that a system based on structural heterogeneity might have the advantage of robustness. Graded release in a stochastic model must depend on a fine kinetic balance between opening of the RyR and, in the presence of sustained Ca^{2+} influx into the subspace, inactivation of the RyRs. As a consequence, conditions such as ischemia, which alters RyR gating (6), could have a disproportionate effect on excitation-contraction coupling. A stochastic system may also be strongly sensitive to the exact nature of single channel kinetics of the L-type Ca^{2+} channel (e.g., burst dwell times), although this is probably less well constrained. In contrast, the heterogeneity model only requires the RyR to have an all-or-none threshold to produce graded release. The structural heterogeneity model indicates that the RyR and L-type Ca^{2+} channel might offer safer targets for pharmacological modification than is suggested by stochastic models. Similarly, Ca^{2+} loading of the SR decreases the threshold subspace Ca^{2+} required to trigger RyR firing, due to an increased Ca^{2+} flux through low probability but non-zero RyR open states. In both models it is easy to see how Ca^{2+} loading in the SR and a slight rise in the myoplasmic Ca^{2+} might give rise to a spontaneous release of Ca^{2+}. In the case of spatial heterogeneity, only a small portion of the subunits (i.e., those with a large number of RyRs) would fire and partially relieve the SR load, leaving the remaining RyRs ready to fire normally.

Ca^{2+} release by the RyR is highly regulated. Perhaps the most important advantage that heterogeneity provides for cell Ca^{2+} regulation is an additional mechanism for regulation of release. For example, Ca^{2+} stored in the lumen of the SR tends to increase the open probability of the RyRs in a cluster (24). The results of our simulations suggest that spatial colocalization may contribute to this phenomenon through increased subspace Ca^{2+} accumulation in subsystems with large numbers of RyRs. Similarly, coupled gating, in which there is a tendency of RyRs to open and close together, may be physical cooperativity, potentially through FK-binding proteins (24). This cooperativity of RyRs in producing sparks may also be augmented by the all-or-none behavior of subunits of a particular size in our model.

Spatial heterogeneity changes in pathological conditions. Hypertrophic heart failure is associated with considerable changes in cell morphology and Ca^{2+} handling (22, 23). Changes in the number and regional localization of L-type Ca^{2+} channels and RyRs also occur. Cell swelling is associated with transient ischemia and may change the physical arrangement of junctional spaces. The reduction in spatial organization may be partially compensated by hyperphosphorylation of the RyR, resulting in a much lower threshold for RyR activation by Ca^{2+} (41). Such a change may compensate for the increased diffusional distances and decreased L-type Ca^{2+} channel fluxes and help maintain the Ca^{2+} release function.

By emphasizing the potential role of spatial heterogeneity in organizing and controlling graded Ca^{2+} release, we do not intend to rule out stochastic mechanisms. They most likely play an important role. Most channel models are at some level stochastic, which reflects the role of Brownian motion in gating. The random events of gating may contribute to the overall graded release response. The lack of some of these random events may contribute to one of the present models inherent limitations. The model, in its current form, does not produce the experimentally observed voltage-dependent “gain” in Ca^{2+} release, i.e., for the same level of Ca^{2+} entry through L-type Ca^{2+} channels, there is more CICR at negative potentials than there is at positive levels. The reason for this phenomenon is that Ca^{2+} release occurs at a relatively distinct threshold. At positive potentials, where the Ca^{2+} channel open probability is high, there is close correlation between the macroscopic current and the maximum current into each subspace. In contrast, at negative potentials, the Ca^{2+} channel open probability is low. The net macroscopic current is a poor approximation of the maximum current influx into a single subunit. Discrete events make Ca^{2+} influx highly variable, and at some point this variation or “noise” has a high probability of achieving the threshold for activation of the Ca^{2+} release unit. Such noise is highest for large unitary currents and low probabilities (62). Therefore, the low voltage gain for Ca^{2+} release could potentially be modeled in this system by weighting the Ca^{2+}-dependent coefficient for opening of the Ca^{2+} release channel by the predicted variance for each population of release subunits. We have chosen not to include this in our model to demonstrate how many of the properties of the Ca^{2+} handling system can be accounted for simply by taking into account anatomical variability and other structural considerations.

We have demonstrated a mechanistically based model of coupling between L-type Ca^{2+} channels and Ca^{2+} release. This model produces graded Ca^{2+} release. It reproduces biexponential Ca^{2+}-dependent inactivation and realistic rates of recovery from inactivation. It does this using experimentally observed variability in microanatomy and a L-type Ca^{2+} channel gating model based on the emerging understanding of structure and function of channel gating in cloned voltage-gated channels. This model is computationally simpler than models requiring Monte-Carlo simulations and extensive supercomputing facilities. Finally, there are testable differences between this and other models concerning the mechanism of graded Ca^{2+} release.

## APPENDIX A: DIFFERENTIAL EQUATIONS

(A1) (A2) (A3) (A4) (A5) (A6) (A7) (A8) (A9) (A10) (A11) (A12) (A13) (A14) (A15) (A16) (A17) (A18) (A19) where *k* = 1, 2,..., *N* (A20) (A21) (A22) (A23) (A24) (A25)

## APPENDIX B: CALCIUM CURRENTS

### L-Type Ca^{2+} Current

(B1) where *N* is the number of calcium release subsystems (B2)

### Ca^{2+} Pump Current

(B3)

### Na^{+}/Ca^{2+} Exchanger Current

(B4)

### Ca^{2+} Background Current

(B5) (B6)

## APPENDIX C: CALCIUM FLUXES

(C1) where *x* is rel, tr, or xfer for release, transfer from the junctional SR to the network SR, and transfer from the subspace to the myoplasm fluxes, respectively.

The fluxes *J _{x}*

_{,}

*are (C2) (C3) (C4)*

_{k}Other fluxes are (C5) (C6) (C7)

## Acknowledgments

GRANTS

This work was supported in part by National Heart, Lung, and Blood Institute Grants R01 HL-59526-01, American Heart Association Established Investigator Award 9940185N, National Science Foundation Knowledge and Distributed Intelligence Grant DBI-9873173, a Whitaker Foundation Biomedical Engineering Research Grant, and a grant from the Oishei Foundation.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2004 by the American Physiological Society