This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.047449v1 87/5/3351    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal 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 Google Scholar Google Scholar Articles by Shannon, T. R. Articles by Bers, D. M. Search for Related Content PubMed PubMed Citation Articles by Shannon, T. R. Articles by Bers, D. M.

A Mathematical Treatment of Integrated Ca Dynamics within the Ventricular Myocyte

Thomas R. Shannon ^* , Fei Wang , José Puglisi , Christopher Weber  and Donald M. Bers 

^* Department of Molecular Biophysics and Physiology, Rush University, Chicago, Illinois; and Department of Physiology, Loyola University-Chicago, Maywood, Illinois

Correspondence: Address reprint requests to Donald M. Bers, E-mail: dbers{at}lumc.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
We have developed a detailed mathematical model for Ca²⁺ handling and ionic currents in the rabbit ventricular myocyte. The objective was to develop a model that: 1), accurately reflects Ca-dependent Ca release; 2), uses realistic parameters, particularly those that concern Ca transport from the cytosol; 3), comes to steady state; 4), simulates basic excitation-contraction coupling phenomena; and 5), runs on a normal desktop computer. The model includes the following novel features: 1), the addition of a subsarcolemmal compartment to the other two commonly formulated cytosolic compartments (junctional and bulk) because ion channels in the membrane sense ion concentrations that differ from bulk; 2), the use of realistic cytosolic Ca buffering parameters; 3), a reversible sarcoplasmic reticulum (SR) Ca pump; 4), a scheme for Na-Ca exchange transport that is [Na]_i dependent and allosterically regulated by [Ca]_i; and 5), a practical model of SR Ca release including both inactivation/adaptation and SR Ca load dependence. The data describe normal electrical activity and Ca handling characteristics of the cardiac myocyte and the SR Ca load dependence of these processes. The model includes a realistic balance of Ca removal mechanisms (e.g., SR Ca pump versus Na-Ca exchange), and the phenomena of rest decay and frequency-dependent inotropy. A particular emphasis is placed upon reproducing the nonlinear dependence of gain and fractional SR Ca release upon SR Ca load. We conclude that this model is more robust than many previously existing models and reproduces many experimental results using parameters based largely on experimental measurements in myocytes.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

Where reliable experimental data exist, modeling data should be as consistent as possible with it. Models will never be as useful as laboratory data in terms of providing concrete information about physiological processes. The complexity of integrated physiological systems means a model will virtually never be perfect. Indeed, it is in this imperfection that we find perhaps a model's greatest usefulness. Where a model cannot be brought to describe experimental data, it identifies areas where our quantitative understanding is incomplete. The identification of these areas coupled with novel mathematical formulation allows the integrated model to be used as a powerful experimental tool to guide investigation and analytical thought.

There are many recent examples of mathematically modeled electrical and ionic homeostasis in the intact cardiac myocyte (McAllister et al., 1975; Beeler and Reuter, 1977; DiFrancesco and Noble, 1985; Priebe and Beuckelmann, 1998; Nordin, 1993; Pandit et al., 2003; Shiferaw et al., 2003; Bondarenko et al., 2004). Among the most valuable of these is the elegant model of Luo and Rudy (1994a,b), which does an outstanding job of modeling the individual currents that combine to produce a credible action potential (AP), and those of the Winslow group (Jafri et al., 1998; Winslow et al., 1999) that were a major step forward in modeling whole-cell Ca homeostasis. Each of these models has its distinct advantages and has proven useful in investigating a variety of different conditions and diseases such as heart failure (Winslow et al., 1999; Puglisi and Bers, 2001).

Since the introduction of many of these models in their original form, new information regarding the processes involved in ECC has become available. For instance, it is now well established that sarcoplasmic reticulum (SR) Ca content has a profound effect upon the process of Ca-induced Ca release (Spencer and Berlin, 1995, 1997; Janczewski et al., 1995; Bassani, J. W. M. et al., 1995; Donoso et al., 1995; Herrmann-Frank and Lehmann-Horn, 1996; Lukyanenko et al., 1996, 2001, 1999; Dettbarn and Palade, 1997; Satoh et al., 1997; Song et al., 1997; Eisner et al., 1998; Hüser et al., 1998; Kurebayashi and Ogawa, 1998; Shannon et al., 2000a, 2002b; Terentyev et al., 2002). This release results in an increase in Ca that appears to be higher just under the sarcolemmal membrane (SL) than in the bulk cytosol and this higher [Ca] is detected by the Na-Ca exchanger and the Ca-dependent Cl channel, which are allosterically activated by it (Hilgemann et al., 1992; Matsuoka et al., 1995; Trafford et al., 1995; Weber et al., 2001, 2002).

With the above advances in mind, we have developed a new model of rabbit cardiac myocyte Ca and Na homeostasis. Our goals were fivefold: i), to track ion influx and efflux such that the model comes to steady state with a realistic balance of Ca fluxes; ii), to incorporate a form of Ca-induced Ca release (CICR); iii), to use realistic parameters that are consistent with laboratory observations; iv), to be able to simulate basic physiological phenomena when all of the components are combined; and v), to incorporate reasonable compromises to allow the model to be solved numerically on a desktop computer.

The model includes the following novel features: i), the addition of a subsarcolemmal compartment along with the junctional and bulk cytosolic compartments to allow proteins in the membrane to sense ion concentrations that differ from bulk; ii), the use of realistic cytosolic Ca buffering parameters; iii), a reversible SR Ca pump as suggested by the results of Shannon et al. (2000a); iv), a [Na]_i-dependent Na-Ca exchanger that is physiologically regulated by Ca as proposed by Hilgemann et al. (1992) (see also Weber et al., 2001); and v), a model of SR Ca release (Fabiato, 1985) including both inactivation/adaptation (Cheng et al., 1995; Stern et al., 1999) and SR Ca load dependence (Shannon et al., 2000b). The model data describe Ca handling characteristics of the cardiac myocyte and the SR Ca load dependence of these processes is accounted for. The model includes a realistic balance of cellular Ca removal mechanisms (Bassani et al., 1994; Puglisi et al., 1999) and the phenomena of rest decay and frequency-dependent inotropy. A particular emphasis is placed upon reproducing the nonlinear dependence of gain and fractional SR Ca release upon SR Ca load (Shannon et al., 2000a). We propose that this model is a robust representation that has useful characteristics not existent in previous models of excitation-contraction coupling (ECC).

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
General approach
All programming was done on IBM PC desktop computers in the standard C programming language. Differential equations were solved using the CVODE package (http://www.netlib.org/ode/index.html). The described system of differential equations will typically generate a Ca transient (stimulation rate 1 Hz) in 2 s on a Dell Precision 3500 workstation with Pentium 4 computer chip with 2.4-GHz processing speed and with 512 MB memory. The computation time was slowed by writing large numbers of parameters to the output text file. All simulations were done running Linux (Red Hat 9.0, Red Hat, Raleigh, NC).

The model is composed of a series of differential equations describing changes in [Ca], [Na], and membrane voltage over time. The AP is generally reconstructed from individual equations representing sarcolemmal membrane channels (Fig. 1) as in Luo and Rudy (1994a,b) with variations in and additions to individual equations and parameters as indicated below (see also Puglisi and Bers, 2001). Markovian channel models were generally avoided unless they were determined to be necessary. Such models often add additional differential equations to the model that must be solved, thus slowing the computation time.

View larger version (37K): [in this window] [in a new window]   FIGURE 1  (A) Diagram of the cell with Ca- and Na-dependent components of the model. (B) Diagrams of the arrangement of junctions over the area of the cell membrane (left) and of the cytosolic space of the sarcomere. Diffusion distances were calculated based upon this arrangement from the center of each compartment assuming that each junction contributed an equal amount of Ca to each space. (C) (Left) An action potential and the accompanying bulk cytosolic Ca transient. (Right) Relevant K and Cl currents. (see Materials and Methods above; Puglisi and Bers, 2001)

View larger version (24K): [in this window] [in a new window]   FIGURE 4  (A) Markovian state model of the ryanodine receptor. The model was developed from that of Stern et al. (1999) with the addition of SR Ca load dependence to the binding of Ca to the activation and inactivation cytosolic sites. (B) SR Ca release flux with RyR-dependent SR Ca leak (inset). Total leak is equal to this diastolic release plus the passive leak flux (4 µM/s). (C) Profile of the four channel states over the course of an AP. (D) Time-dependent profile of the channel open state with changes in the kon for Ca binding with [Ca]SR (inset).

This spacing of junctions also gives a distance of diffusion to the middle of the SL compartment of 500 nm. This is a narrow space just under the SL membrane that is diffusionally restricted on the cytosolic side by such structures as the mitochondria and the myofilament array. Similarly, an average distance of diffusion of 0.45 µm from the SL to the middle of the half sarcomere (i.e., the cytosolic compartment) was also calculated assuming a sarcomere length of 1.8 µm. These diffusional distances dictate the transfer functions for the fluxes between each compartment. This adjusts the lumped nature of the compartments to an appropriate physical geometry.

Na buffers were located only in the junction and SL compartment and were modeled as rapidly binding molecules with the standard Hill equation (Eq. 14, Table 3; parameters were taken from Bers et al., 1986).

(1)

(2)

(3)

(4)

(5)

(6)

where I_Na is fast Na current (Eq. 37), I_NaBk is Na leak (Eq. 38), I_NaK is Na-K pump current (Eq. 41), I_Ca is L-type Ca current (Eq. 93), I_CaBk is SL Ca leak (Eq. 94), I_NCX is Na-Ca exchange current (Eq 96), and I_SLCaP is SL Ca pump transport (Eq. 97).

I_NaBk and I_NaK were formulated as previously (Luo and Rudy, 1994a,b; Puglisi and Bers, 2001) and were adjusted to give resting flux values similar to those measured experimentally by Despa et al. (2002a,b).

Ca currents
Ca current was formulated as previously described (Luo and Rudy, 1994a,b; Puglisi and Bers, 2001) with modifications. The formulations are based upon the Goldman-Hodgkin-Katz equation. The permeabilities of the L-type Ca channels for Na and K were decreased, consistent with the selectivity of the channel for these ions in the presence of Ca in the perfusate (Hess et al., 1986; Tsien et al., 1987). Currents were scaled to approximate those of Puglisi et al. (1999). Ca-dependent inactivation was modified to be calmodulin dependent (Eq. 77; Peterson et al., 1999; Qin et al., 1999). Kinetics of binding were set with a K_d of 7 µM to give typical experimental values for the fast and slow time constants of decline in I_Ca due to inactivation during a normal Ca release (3 and 28 ms, respectively), when a depolarization from –40 to 0 mV was simulated (Eq. 78). Ninety percent of the channels are located in the junctional cleft membrane (Scriven et al., 2000). T-type Ca current was not included in the model because this current, although present in some cardiomyocytes (Bean, 1989; Nuss and Houser, 1993), is undetectable in normal rabbit ventricular myocytes (Yuan et al., 1996).

The SL Ca leak flux was adjusted to match that calculated from Negretti et al. (1993). The SL Ca pump was formulated using the standard Hill equation and approximates cellular estimates of Bassani et al. (1994).

The Na-Ca exchanger was formulated essentially as described in Weber et al. (2001). The scheme is a general improvement over other models in that it is [Na]_i dependent and allosterically regulated by Ca. As with all of the channels and transporters with the exception of the L-type Ca channel, Na-Ca exchange molecules were evenly distributed throughout the cell membrane, 89% in the SL compartment and 11% in the junctional membrane.

Sarcoplasmic reticulum Ca transport
Ca transport out of the SR into the cytosol was calculated:

(7)

where J_rel is the SR Ca release flux (Eq. 106) and J_pump is the SR Ca pump flux (Eq. 98). The SR Ca pump is a reversible enzyme formulated as in Shannon et al. (2000a). It is important to note that the net SR Ca pump rate never reverses. The enzyme merely slows as product ([Ca]_SR) builds up. All of the SR Ca pump molecules were localized to the cytosolic compartment (Fig. 1).

The cytosolic Ca dependence of the SR Ca release channel (RyR) was modeled in a steady-state manner as in Stern et al. (1999) (see Fig. 4 A, based upon Fabiato, 1985 and Cheng et al., 1995). The model has four states: resting or closed (R), open (O), inactivated (I), and resting inactivated (RI). Models of this type are typically used by investigators working on single-channel RyR gating (Scheifer et al., 1995; Zahradnikova and Zahradnik, 1996; Keizer and Levine, 1996; Villalba-Galea et al., 1998). Additional states are often added to such models but this one with our modifications (e.g., lumenal [Ca]_SR dependence) works well.

During a typical Ca release cycle, the RyR exists predominantly in the closed state. Upon depolarization and sudden Ca influx into the cleft, Ca binds to the cytosolic activation site on the RyR resulting in a transition to the open state. Ca also binds to the slower cytosolic inactivation sites resulting in the transition from the open state to the inactivated state. As Ca declines, there is a transition to the resting inactivated state and finally back to the closed state.

SR Ca load dependence was added as a lumenal Ca binding site that affects the state of the RyR. Ca binding to this site affects both the inactivation and activation cytosolic Ca-dependent rate binding constants (Eqs. 99–101). Neither rate constant was allowed to exceed diffusional limitations for Ca binding. When Ca is bound on the lumenal side and this site is occupied, the affinity of the cytosolic activation site increases whereas the affinity of the inactivation site decreases. Thus, there is a greater tendency for the RyR to undergo transition to the open state. On the other hand, when less Ca is bound to the lumenal site, there is a greater tendency for the RyR to exist in the inactivated state.

All of the SR Ca release takes place in the junctional compartment. A passive leak from the SR into the junctional compartment was also added to give a total diastolic leak (passive and diastolic RyR leak) of 4 µM/s (Shannon et al., 2001).

Other currents
Most K currents were formulated as by Luo and Rudy (1994a,b) as modified by Puglisi and Bers (2001) and Bassani et al. (2004). Two notable exceptions are in the I_Ks current and the Ca-dependent Cl current (I_Cl). I_Ks was modeled to resemble the data of Tohse (1990). In accordance with this article, currents were made Ca dependent and were scaled as appropriate for the rabbit. Cl currents were added to the model and scaled to approximate those of Zygmunt and Gibbons (1991) and Puglisi et al. (1996) using a K_d for Ca of 100 µM (Collier et al., 1996)

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
During excitation-contraction coupling, Ca enters through Ca channels in response to an action potential (Figs. 1 and 2). The AP in our model is formulated by combining equations that describe the characteristics of the voltage-dependent ion channels in isolation into a coherent model. The equations used for the Na and K currents are based upon those of Luo and Rudy (1994a,b) with some modifications (Jafri et al., 1998; Puglisi and Bers, 2001; Bassani et al., 2004). Combined, the currents produce a relatively normal AP waveform with a duration of 235 ms.

View larger version (31K): [in this window] [in a new window]   FIGURE 2  (A) Diagram of the way in which the L-type Ca channel operates. The model includes Ca-dependent inactivation that is Ca-calmodulin (CaM) dependent. Most of the channels (90%) are in the cleft space (Scriven et al., 2002). (B) L-type Ca current generated during the AP (right). The currents of Puglisi et al. (1999) are shown for comparison (left).

In our model 90% of the Ca channels are located in the junction based upon the data of Scriven et al. (2000, 2002) who found that the vast majority of the channels colocalize with the RyR. AP-dependent Ca currents tend to be of variable shape (Fig. 2 B; Linz and Meyer, 1998; Puglisi et al., 1999; Linz and Meyer, 2000). This variability is species dependent and may be due to changes in the shape of the AP, to the effects of temperature, and to the effects of released Ca upon the current. Ca-dependent inactivation of the L-type Ca channel is accounted for within the model and changes in the released Ca also cause the expected changes in the shape of the Ca current (Fig. 5 B; see below). The shape of the current in the model with a normal release during 1-Hz pacing is typical (Fig. 2 B; see also Fig. 5 B).

View larger version (25K): [in this window] [in a new window]   FIGURE 5  (A) CICR in response to graded L-type Ca current. The change in the bulk Ca during the transients and the integrated SR Ca release (top) correspond to the peak L-type Ca currents (bottom) at each voltage. (B) L-type Ca current in response to a 200-ms square pulse to 0 mV from a holding potential of –80 mV. Currents show normal Ca-dependent inactivation of current in the presence of release and with Ca buffered to 100 nM when compared to "EGTA" currents with no Ca-dependent inactivation. ICa and integrated SR Ca release are normalized to their peaks (0.410 µM and 57.4 µmol/l cytosol, respectively). (C) [Ca]SRT dependence of the gain (left) and fractional release (right). The data of Shannon et al. (2000a) are shown for comparison (•).

View larger version (26K): [in this window] [in a new window]   FIGURE 3  (A) Ca buffers in each compartment. (B) Bulk cytosolic Ca transient. (C) Ca transients in each of the three non-SR compartments. Data are compared to that of Weber et al. (2002).

The Ca dependence of release was tested by simulating a standard experimental assay for ECC. The digital cell was stimulated to steady state with 200-ms square pulses to 0 mV and a test pulse was applied to varying voltages from –60 to +40 mV. The release of Ca from the SR was graded and dependent upon the L-type Ca current (Fig. 5 A).

Importantly, released Ca also has a major effect upon the L-type Ca current. As would be expected in an isolated myocyte, this current exhibits Ca-dependent inactivation in the model. This inactivation is known to be dependent upon Ca binding to calmodulin (Peterson et al., 1999; Qin et al., 1999). Fig. 5 B graphically shows the effect of this Ca binding in the model. Ca current with normal release into the cleft inactivates rapidly as would be expected. On the other hand, Ca current with no release and Ca current with [Ca]_jct clamped at 100 nM inactivates slower due to reduced Ca binding to the channel/calmodulin complex.

Similarly, Ca release was tested as a function of [Ca]_SRT. Experiments by Shannon et al. (2000a) were simulated. After depletion of SR Ca, the cell was depolarized to 0 mV for 200 ms per pulse in the absence of Na inside and outside of the cell. The absence of Na prevented Na-Ca exchange from operating, which causes the retention of nearly all Ca brought in by the L-type Ca current within the cell. Most of the Ca is retained within the SR resulting in a larger release with each successive depolarization. The model simulates the actual data indicating normal [Ca]_SRT-dependent release (Fig. 5, C and D).

The overall physiology of the release from the SR is also in agreement with laboratory data. J. Bassani et al. (1993b, 1995) found that the fractional release from the SR during the [Ca]_i transient is 50%, which agrees well with the model as do the [Ca]_SR transients of Shannon et al. (2003a) (Fig. 6, A and B). Fig. 6 A shows the effect of removing the lumenal regulation from the SR Ca release process. The first beat upon removal of the regulation has a larger fractional release with a longer time to peak.

View larger version (24K): [in this window] [in a new window]   FIGURE 6  (A) Changes in [Ca]SRT and [Ca]SR over the time course of the Ca transient. The fractional release (30%) agrees well with the value measured by J. Bassani et al. (1993b, 1995). The [Ca]SR transient agrees well with the data of Shannon et al. (2003) (B) Net SR Ca pump flux (left) into the SR, which is the difference between the forward and reverse unidirectional fluxes (right) (Shannon et al., 2003a).

The formulation of the Na-Ca exchanger is based upon that used in Weber et al. (2001) and the modeled current during the AP is similar to the actual data (Fig. 7). In contrast to previous formulations used in whole-cell simulations (Luo and Rudy, 1994a; Jafri et al., 1998), the equation accounts for the [Ca]_SL-dependent activation of Na-Ca exchange that takes place during the Ca transient. Note once again that this activation depends upon [Ca]_SL, i.e., the [Ca] just under the SL membrane that is actually sensed by the exchange protein. The equation also accounts for the effects of intracellular [Na] at the transport site, specifically Na just underneath the SL membrane in the SL compartment, another important addition.

View larger version (20K): [in this window] [in a new window]   FIGURE 7  (A) Diagram illustrating the operation of Na-Ca exchange within the model. The exchanger includes an allosteric Ca-dependent activation (Hilgemann et al., 1992; Weber et al., 2001). (B) Model exchange current during an AP (bottom). The data of Weber et al. (2002) are shown for comparison (top).

The transport competition among Na-Ca exchange, the SR Ca pump, and the SL Ca pump takes place in a characteristic, species-dependent way (Bassani et al., 1994). The balance of integrated fluxes in rabbit cardiac myocytes is depicted in Fig. 8 and the model simulates this data quite well. The data indicate that the transport of Ca from the SL and bulk cytosolic compartments is correctly balanced and therefore imply that the balance during ECC and release at steady state is also correct. Note that these integrated steady-state Ca fluxes during a twitch also define the gain of ECC (integrated Ca released/integrated Ca influx) in rabbit during an action potential as 3–4, consistent with cellular data.

View larger version (16K): [in this window] [in a new window]   FIGURE 8  The balance of fluxes that transport Ca from the cytosol. Each transporter competes for Ca resulting in a defined percentage of the total Ca translocated by each. The data of Puglisi et al. (1996) are provided for comparison (left). The of [Ca] decline for the model experiments (190 ms normal twitch, 750 caffeine alone, and 7000 for 0 Na/0 Ca plus caffeine) were similar to the experimentally obtained values there.

View larger version (29K): [in this window] [in a new window]   FIGURE 9  Response to a sudden increase (red) or decrease (blue) in [Ca]SRT (bottom). The model reacts by increasing the release from the SR, thus increasing the peak of the bulk Ca transient (top) and decreasing the [Ca]SRT for the next beat until the SR Ca content returns to normal and the model returns to the previous steady state.

The physiological utility of the model was tested by simulating well-known cellular responses. The first of these was the phenomenon known as rest decay (Fig. 11 A; Bers et al., 1993). The cell was stimulated to steady state and then stimulations were suddenly terminated. Upon resting the cell, the SR lost Ca at a rate comparable with that seen in normal rabbit myocytes. Moreover the decline in amplitude of postrest twitch Ca transients also matched rest-decay data (Fig. 11 B). The larger percent decline in twitch [Ca]_i as a function of [Ca]_SRT is a consequence, at least in part, of the steep lumenal Ca dependence of the RyR gating (Fig. 5). Note well that because this relationship levels off before it reaches 0 SR Ca, the SR does not completely drain. Similar to the resting rabbit ventricular myocyte, the resting digital cell has Ca available for release.

View larger version (21K): [in this window] [in a new window]   FIGURE 11  The model exhibits normal rest decay (A). Summary data for both [Ca]SRT and twitch [Ca]i from the model (bottom) are to be compared with the data of Bers et al. (1993) (top).

View larger version (26K): [in this window] [in a new window]   FIGURE 10  The model exhibits a positive force-frequency relationship in response to an increase in frequency from 0.125 to 0.5 Hz (A). Panel B shows summary data from the model (bottom) along with the corresponding rise in diastolic Na concentration in the bulk cytosol (inset). Data from Maier et al. (2000) are shown for comparison.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
This article describes a new model that describes the Ca dynamics within the cardiac myocyte. The model was designed to be simple yet descriptive enough to simulate complex physiological responses.

CICR
The goals of the project included the use of CICR to initiate and control SR Ca release while keeping the model available for use on a personal desktop computer with reasonable computational power. This was done by avoiding computationally expensive formulations that include extensive, spatially descriptive, local junctional, control and/or Monte Carlo simulations of stochastic RyR properties. The model uses a "bulk" junctional compartment that is a small but realistic proportion of the total cell volume (<0.1%). The gating of the SR Ca release channels is graded (Fig. 5) and dependent entirely upon the [Ca] within this cleft (with important modulation of that dependence by [Ca]_SR). This makes the entire process of ECC more realistic and enables more effective investigation of phenomena that affect CICR.

Although the simplified control via CICR is quite good, the description in this model is still limited. For instance, SR Ca release as a function of voltage has been shown in the elegant experiments of Sipido, Balke, and Weir (Sipido and Wier, 1991; Balke et al., 1994) to have a peak that is slightly shifted to the left of that of the bulk L-type Ca current. The model as it is currently formulated does not include this characteristic. Neither this model, nor any model that lacks the necessary spacial resolution within the cleft is likely to describe such a characteristic accurately.

The RyR is modeled with a modification of the scheme best described by Stern et al. (1999) (Fig. 4). The scheme works remarkably well as described with few changes in the parameters they used. The peak open probability of 1% is reasonable when data from the literature are considered. There are 1.5 million RyR per myocyte (Bers and Stiffel, 1993; Bers, 2001). For a single RyR flux of 0.4 pA (Majia-Alvarez et al., 1999) and a peak SR Ca release flux of 2 mM/s, only 2% of the RyR would be open.

Importantly, however, we have added dependence of those parameters upon [Ca]_SR (Fig. 4). Regulation of SR Ca release by SR Ca is well documented (Bassani, J. W. M. et al., 1995a; Janczewski et al., 1995; Eisner et al., 1998; Lukyanenko et al., 1999; Shannon et al., 2002) and may be mediated by a lumenal Ca sensor (Györke and Györke, 1998; Györke et al., 2004).

As formulated, the [Ca]_SR affects the RyR in a classic Michaelis-type relationship (Eq. 99). Lower [Ca]_SR does not directly inactivate the RyR. Instead, lumenal Ca alters the Ca binding parameters on the cytosolic side of the membrane such that both the lumenal and junctional [Ca] act in combination to both initiate and terminate SR Ca release. The model accounts for the dependence of the ECC process upon SR Ca load largely because of this modification (Fig. 5). Notably it explains the failure of CICR when [Ca]_SR is <40–50% of its normal value (Bassani, J. W. M. et al., 1995a; Shannon et al., 2000). It also explains experimental data that have led to conclusions that lumenal Ca depletion plays a key role in the termination of release, not by exhausting the Ca gradient per se, but rather by an allosteric regulatory effect. This extends the list of possible mechanisms for interruption of positive feedback of CICR that were considered by Stern et al. (1999).

Ca transporters
Four Ca transport systems can compete for Ca during Ca transient decline: i), SR Ca ATPase; ii), SL Na-Ca exchange; iii), SL Ca ATPase; and iv), mitochondrial Ca uptake. The first three are well described in the model whereas the last (which is minor) will be added at a future date.

The Na-Ca exchanger undergoes activation upon Ca binding to an allosteric site located on the exchanger on the intracellular side of the SL membrane (Hilgemann et al., 1992; Matsuoka et al., 1995; Trafford et al., 1995; Weber et al., 2001, 2002). This site has been shown to be physiologically relevant with measurements of the K_d for Ca varying from 125 to 450 nM (Weber et al., 2001).

The new formulation of the Na-Ca exchanger is based upon that in Weber et al. (2001, 2002), which matched extensive data therein. Importantly, the equation accounts for the deactivation of Na-Ca exchange that takes place during the Ca transient as [Ca]_i (and [Ca]_SL) falls to rest and Ca dissociates from its allosteric activation site. The addition is a double-edged sword. It allows the Na-Ca exchanger to turn off near resting [Ca]_i, preventing the high rate of Ca loss that would otherwise continue, making diastolic [Ca]_i too low, and depleting the cell of Ca, yet allows for rest decay (Fig. 10; Bers, 1985; Dani et al., 1979; McCall et al., 1996; Satoh et al., 1997; Terracciano et al., 1995). However, this allosteric Ca-dependent activation is also responsible for a higher than experimentally assessed initial peak in inward Na-Ca exchange current (Fig. 7 B; Weber et al., 2001) that corresponds to the peak of the Ca transient in the junction and along the SL. Ongoing investigation into the time dependence of this activation may yield valuable results that may explain this apparent dichotomy and allow a more accurate simulation of the Na-Ca exchange current in the future.

Additionally, unlike some previous formulations (Luo and Rudy, 1994a), the equation also accounts for the effects of intracellular [Na], which is often reduced or removed experimentally to alter the exchange activity. In addition, [Na]_i is known to change within the myocyte during changes in frequency and under the condition of chronic heart failure (Despa et al., 2002b), and this surely is a factor contributing to the force-frequency relationship as shown in the inset of the bottom panel of Fig. 10 B. Na also rises in the other compartments to a similar extent. Future studies both in the model as well as at the bench should provide useful quantitative information about the behavior of this important ion.

The SR Ca pump has also been reformulated in a more realistic way. This ATPase is traditionally represented in a unidirectional Hill-type relationship (Sipido and Wier, 1991; Luo and Rudy, 1994a,b; Balke et al., 1994; Bassani et al., 1994). However, it is well known that the pump is a reversible enzyme being capable of mediating both forward and reverse unidirectional fluxes, although the net rate never reverses physiologically (Weber et al., 1966; Makinose, 1971; Takenaka et al., 1982; Feher, 1984; Shannon et al., 2000b, 2001). The extent that the reverse flux takes place to slow uptake under physiological conditions is debatable but by using a formulation that allows for this reversal flux, we have allowed for these potential effects. This representation will also allow for realistic extension of the model into energetics and ATP consumption because net ATP hydrolysis is reduced, when compared to a simpler pump-leak balance.

Critical to the process of SR Ca uptake was the necessity of formulating a correct balance between the SR Ca pump and the leak of Ca from the organelle. The proper quantity of leak was formulated through the combination of RyR flux plus a passive leak for Ca from the SR into the junctional space. The background leak of Ca from the SR of 3 µM/s at normal [Ca]_SRT increased the diastolic value of [Ca]_jct and thus sustained an increase in the diastolic RyR activity that added an additional 1 µM/s of leak to total 4 µM/s SR Ca leak. This value is within the normal range (Fig. 4 B; inset, Györke et al., 1997; Shannon et al., 2001, 2003b). This leak was important in contributing to both the positive force-frequency relationship and rest decay.

Perhaps the most important aspect of the Ca transport is the balance of fluxes (Fig. 8). Bassani et al. (1992) initially described the competition between the relevant transporters for Ca during transient decline. The effectiveness of each of these transporters in this competition varies with species but the human is most like the rabbit, which is described here (Bassani et al., 1993b, 1994; Bassani, R. A. et al., 1995; Negretti et al., 1993; Bers, 2001; Piacentino et al., 2003).

The model is formulated to maintain this balance of fluxes. Of the four relevant transporters only the mitochondrial uptake, a relatively minor flux, is missing. This SL Ca pump transport in the model represents the "slow" transport component typically observed in physiological experiments (Bassani et al., 1993a). The mitochondrial part of this component, 1% of the total Ca efflux from the cytosol during the transient, is distributed among the three existing model fluxes.

The characterization of the fluxes and maintenance of this balance is absolutely critical to the predictive value of the model. For instance, at steady state the amount of Ca that enters the cell is equal to the amount which leaves. Carried further, this means that the amount of Ca taken up into the SR is equal to the amount of release whereas the amount of Ca transported out of the cell on the SL Ca pump and the Na-Ca exchanger is equal to the amount of influx primarily through the L-type Ca channel (plus a minor contribution from outward Na-Ca exchange current). The balance of fluxes indicates the "gain" of ECC as defined by the integrated release divided by the integrated current (i.e., the stimulus of the CICR).

In addition, to extend the simulations into situations where transport of one type is altered, one must predict the extent to which other transporters will compensate. For instance, the balance between these fluxes also varies with pathological conditions such as heart failure (Pogwizd et al., 1999, 2001; O'Rourke et al., 1999; Shannon et al., 2003b), which involves an increase in Na-Ca exchange activity, a decrease in SR Ca pump activity, and an increase in SR Ca leak, depending upon the failure model. The extent to which either or all of these alters [Ca]_SRT depends upon the extent to which the alteration of each affects transport relative to the other and the competition among the transporters must be correct to predict the overall effects (Shannon et al., 2003b).

The subsarcolemmal Ca compartment
Important for the development of both this aspect of the model as well as the currents across the sarcolemmal membrane was the addition of a subsarcolemmal compartment. During the twitch Ca transient there are spatial gradients near the SL membrane. These gradients are likely due to diffusion restrictions that are caused by the mitochondria and the paracrystaline array in the bulk cytosol. Ca-dependent proteins that are located in or near the membrane are therefore likely responding to a [Ca] that is very different from that which is measured in the bulk cytosolic compartment. The presence of the SL compartment is also suggested by the published values for the Ca affinity of membrane-associated Ca-dependent proteins such as the Na-Ca exchanger and the Ca-dependent Cl channel. These are generally lower than would be expected if these proteins were dependent upon the bulk cytosolic Ca, and the SL compartment is needed to appropriately incorporate these proteins into the model. Indeed, it was by measuring Na-Ca exchange activity and I_Cl(Ca) that the difference in the concentration near the SL membrane was inferred by Trafford et al. (1995) and estimated by Weber et al. (2002) (Fig. 3). We have simulated the effects of the higher [Ca] near the SL membrane by creating a new compartment that is separate from that of the bulk cytosol. It is perhaps worth emphasizing that the authors are not actually proposing the existence of a new anatomical compartment within the cell. Rather, this change in the model reflects the existence of a functional compartment where a different Ca profile exists. Although the depth of this compartment is somewhat arbitrary, it was chosen to allow for a physiologically relevant diffusion coefficient for Ca between the SL compartment and the bulk (Table 2). The model accurately reflects the Ca changes within this SL compartment and therefore is an advance over previous ventricular models in that it better accounts for the activity not only of the Na-Ca exchanger but of all Ca-dependent SL proteins. These include the SL Ca pump, the Ca-dependent Cl current, and the L-type Ca channels within this compartment (which undergo Ca-dependent inactivation). The addition of this compartment with a [Ca] that is generally higher than the bulk may be important for accurate simulation of such phenomena as delayed after depolarizations, which are often related to Na-Ca exchange-mediated diastolic membrane depolarization (Pogwizd et al., 2001).

Steady-state behavior
The model described here comes to steady state. This simply means that it continues to operate with repetitive stimulation in a realistic manner such that the cell does not gain or lose Ca over time unless the cellular environment changes. Moreover, the characteristics of the model at steady state are consistent with the normal cellular physiology (e.g., the balance of fluxes and the gain are correct). The model is, therefore, appropriate for use in situations where continued operation for long periods is necessary. This ability is necessary to accurately simulate physiological and pathophysiological reactions, both in the laboratory where many experiments are conditioned with a number of stimulations, and as part of the larger, more integrated context of the whole tissue or organism (Usyk and McCulloch, 2003).

Importantly, the steady state that is achieved is independent of the history of the digital cell. After perturbation of cellular ionic concentrations, the model always returns to the same steady state that is dependent strictly upon the state of the individual components of the ECC process and upon the extracellular environment (Fig. 9). Such behavior allows simulation of physiological responses that are certain to depend strictly upon the experimental perturbation with no unrelated artifacts.

The model accurately reproduces several well-known physiological responses that are made possible by this steady-state behavior. In particular, it reproduces the normal increase in Ca transient peak in response to an increase in frequency (Fig. 10).

In conclusion, we have formulated a realistic mathematical model of Ca homeostasis within the cardiac myocyte. This model uses recent advances in the study of cardiac ECC to improve over previously proposed models. It should be useful for simulating a variety of cellular responses that are demonstrable in the laboratory and will allow more in-depth investigation into the mechanisms behind these responses, particularly when actual experimental data are difficult to obtain.

	APPENDIX

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Environment
All tabulated parameters are for temperature (T) = 310 K (37°C).

Temperature
All parameters and fluxes (J) that were made temperature sensitive within the model were multiplied by a proportionality constant (Q_J) such that:

(8)

where the Q₁₀ is a constant that represents the change in the relevant parameter for each 10°C change in temperature.

Diffusion
Diffusion is considered to take place, on average, from the center of each compartment to the center of its adjacent compartment. The junctional compartment communicates only with the SL compartment, the SL compartment with the junction and the bulk cytosolic compartments, and the bulk cytosolic compartment only with the SL.

For ion I diffusion over a distance x through area A from the center of compartment C1 to the center of compartment C2 (Table 2):

(9)

Nernst potential
For each compartment C, the Nerst potential E_C is calculated as:

(10)

(11)

(12)

(13)

Voltage
Change in voltage was considered to be the result of the sum of all currents across the membrane and no regional variations in voltage were considered.

Na buffering
Significant, rapid Na buffering was only considered to be associated with cell membrane and therefore was only present in the SL and junctional compartments. The change in Na binding to ligand L over time (d[Na·L]/dt) for each compartment was calculated as (Table 3):

(14)

Na currents
Na channels and the Na-K pump were all considered to be evenly distributed throughout the cell membrane. The equations are as in Luo and Rudy (1994a, Table 4). The Nernst potential (E_Na) was therefore considered to be the result of the difference between the extracellular compartment and the interior cell compartment under the region of the cell membrane where the channels are present. For instance, junctional Na channels (11% of the total) were dependent upon extracellular [Na] and junctional [Na].

(15)

(16)

(17)

(18)

If V < –40 mV

(19)

(20)

(21)

(22)

(23)

(24)

(25)