Biophysical Journal 85:3666-3686 (2003)
© 2003 The Biophysical Society
Model of Intracellular Calcium Cycling in Ventricular Myocytes
Y. Shiferaw *,
M. A. Watanabe *,
A. Garfinkel 
,
J. N. Weiss and
A. Karma *
* Department of Physics and Center for Interdisciplinary Research on Complex Systems, Northeastern University, Boston, Massachusetts 02115; and Department of Medicine and Cardiology, University of California, Los Angeles, California 90095
Correspondence: Address reprint requests to Alain Karma, Physics Dept., 110 Forsyth St., Northeastern University, Boston, MA 02115. Tel.: 671-373-2929; Fax: 671-373-2943; E-mail: a.karma{at}neu.edu.
 |
ABSTRACT
|
|---|
We present a mathematical model of calcium cycling that takes into account the spatially localized nature of release events that correspond to experimentally observed calcium sparks. This model naturally incorporates graded release by making the rate at which calcium sparks are recruited proportional to the whole cell L-type calcium current, with the total release of calcium from the sarcoplasmic reticulum (SR) being just the sum of local releases. The dynamics of calcium cycling is studied by pacing the model with a clamped action potential waveform. Experimentally observed calcium alternans are obtained at high pacing rates. The results show that the underlying mechanism for this phenomenon is a steep nonlinear dependence of the calcium released from the SR on the diastolic SR calcium concentration (SR load) and/or the diastolic calcium level in the cytosol, where the dependence on diastolic calcium is due to calcium-induced inactivation of the L-type calcium current. In addition, the results reveal that the calcium dynamics can become chaotic even though the voltage pacing is periodic. We reduce the equations of the model to a two-dimensional discrete map that relates the SR and cytosolic concentrations at one beat and the previous beat. From this map, we obtain a condition for the onset of calcium alternans in terms of the slopes of the release-versus-SR load and release-versus-diastolic-calcium curves. From an analysis of this map, we also obtain an understanding of the origin of chaotic dynamics.
|
INTRODUCTION
|
|---|
The contraction of a cardiac myocyte is triggered by an intracellular rise in calcium concentration that is due to a coordinated release of calcium from the sarcoplasmic reticulum (SR) (Fabiato, 1983
). The release of calcium from the SR occurs via ryanodine receptors (RyR), which are in close proximity to L-type calcium channels that are located in the cell surface membrane and T-tubules (Meissner, 1994; Wang et al., 2001). When the cell is depolarized, L-type channels open and allow calcium entry into a confined microdomain. The rise of calcium in this small space is sensed by the nearby cluster of RyR channels that in turn open via calcium-induced calcium release (CICR) (Fabiato, 1983). As the calcium concentration in the cell rises, contractile elements are activated and the cell contracts. An uptake pump, which is activated by the rise in calcium, then pumps calcium back into the SR. This interplay between voltage across the cell membrane and intracellular calcium cycling forms the basis of excitation-contraction (EC) coupling.
During normal beating of the heart, myocardial cells undergo periodic depolarizations of the membrane called action potentials (AP). The shape of the AP waveform is determined by the flux of ions across the membrane. Some of these fluxes, such as those due to the L-type channel current (ICa) and the NaCa exchange current (INaCa), are modulated by intracellular calcium concentration. Thus, the calcium system is driven by an AP waveform that is itself dependent on the dynamics of the calcium system. Recently, Chudin et al. (1999) were able to shed light on this coupling by stimulating a rabbit ventricular myocyte using a clamped AP waveform. They observed that at high stimulation rates the whole cell calcium transients exhibited alternans even when the clamped voltage stimulus was periodic. This result demonstrates that the calcium system could be dynamically unstable independently of the dynamics of the membrane potential. Thus, given the bidirectional coupling between voltage and calcium, this raised the important possibility that abnormalities in calcium cycling could influence the membrane potential and hence promote arrhythmias. For example, ventricular fibrillation is typically preceded by ventricular tachycardia, where the ventricles can be driven at rapid rates by an unstable rotating spiral wave. It is possible that at high rates a dynamical instability of calcium cycling at the single cell level may promote wave break and thus underlie the transition of ventricular tachycardia to ventricular fibrillation. Therefore, a physiologically based model of calcium cycling is necessary to investigate this important bidirectional coupling of voltage and calcium. In particular, it is important that this model reproduces the experimentally observed complex dynamics of the calcium system at high stimulation rates.
There have been many efforts to model intracellular calcium dynamics. Initial modeling attempts were faced with the challenge of explaining the experimentally observed (Wier et al., 1994) linear dependence of calcium release from the SR on the whole cell ICa current (graded release), given that RyR channels are activated by a CICR mechanism that favors an all-or-none response. This issue was resolved by Stern (1992), who pointed out that graded release can be explained by the stochastic recruitment of independent local release fluxes. Stern referred to models of this type as local control models since a small cluster of L-type channels were coupled to a cluster of RyR channels via a local pool of calcium. Stern emphasized that models that couple whole cell ICa to RyR flux via a single pool, which he referred to as common-pool models, would not be able to explain graded release. This hypothesis was later confirmed by high resolution confocal imaging, which showed that the global rise in calcium was due to the summation of many local release events sometimes referred to as calcium sparks (Cannell et al., 1994; Cleemann et al., 1998; Lopez-Lopez et al., 1995; Niggli, 1999; Wier and Balke, 1999).
The local nature of the release has led researchers to explore the detailed calcium dynamics within the dyadic junction, i.e., the space between the L-type channel clusters and the corresponding RyR receptor clusters (Cannell and Soeller, 1997; Greenstein and Winslow, 2002; Rice et al., 1999; Peskoff and Langer, 1998; Sobie et al., 2002; Stern et al., 1999). However, it turns out to be difficult to characterize experimentally the gating kinetics of the RyR clusters and the interaction of these clusters with the L-type channels. As a result, these kinetics are still poorly understood and, in particular, the detailed mechanism by which release from the SR is terminated remains controversial (Bers, 2001; DelPrincipe et al., 1999; Niggli, 1999). Moreover, even if the calcium dynamics at the dyadic junction is precisely known, there are on the order of 104 junctions within a cell (Cleemann et al., 1998). Given that each junction is on a submicron scale, a modeling approach that attempts to resolve the dynamics of each junction is unlikely to be computationally tractable for modeling arrhythmogenesis at the tissue or organ level.
Several authors have proposed tractable phenomenological models of calcium cycling at the whole cell level (Chudin et al., 1999; Fox et al., 2002; Glukhovsky et al., 1998; Jafri et al., 1998; Luo and Rudy, 1994; Snyder et al., 2000; Viswanathan et al., 1999). However, several of these models (Jafri et al., 1998; Glukhovsky et al., 1998; Snyder et al., 2000) are common-pool models and, as Stern (1992) predicted, do not reproduce graded release. Hence, they do not correctly describe the coupling between the L-type calcium current and the release from the SR that is crucial to describe EC coupling. Viswanathan et al. (1999) presented an improved version of a model of release flux from the SR first presented in the Luo-Rudy II model (Luo and Rudy, 1994), which incorporates graded release. This model was able to reproduce calcium alternans in simulations performed in a ring of coupled cells, where the voltage waveform was not clamped (Hund et al., 2000). However, this model does not take into account the fact that the release from the SR is the summation of discrete release events. Consequently, it cannot, for example, describe the influence of the amplitude and duration of these events on the dynamics of calcium cycling and the instability mechanism responsible for alternans.
In this article, we present a new model of EC coupling in ventricular myocytes. This model distinguishes itself from previous models in that it represents release from the SR as a summation of elementary release events that correspond to calcium sparks. Furthermore, the model formulation is based on two experimentally measurable constitutive relations:
- The relation between the rate of spark recruitment and the whole cell ICa. This relation describes quantitatively the crucial coupling between ICa and calcium release from the SR. This relation is formulated by introducing a variable, N(t), which corresponds to the total number of sparks recruited at a given time t in the whole cell, and by making the rate of spark recruitment (dN/dt) proportional to the whole cell ICa. This linear relationship between dN/dt and ICa is based on the available experimental data on simultaneous measurements of spark occurrence and whole cell ICa (Collier et al., 1999; Lopez-Lopez et al., 1995) and reproduces graded release.
- The relation between the total amount of calcium released form the SR and the SR load. This relation can be built into the model by letting the rate of spark recruitment and/or the average release current of a single spark depend on the junctional SR (JSR) concentration. This dependency is then chosen phenomenologically to reproduce the highly nonlinear relationship between fractional release and load that has been measured experimentally (Bassani et al., 1995).
Another important experimental input into the model is the observation (Cleemann et al., 1998; Niggli, 1999) that the lifetime of sparks is relatively constant and independent of both SR load and diastolic calcium level. Based on this observation, we model termination of release from the SR by assuming that the release current of a single spark decays exponentially in time from its initial peak value, with a time constant that is taken comparable to the observed spark lifetime.
The model includes additional features that are not present in previous models. The first is a phenomenological description of intracellular sodium accumulation. This accumulation is well documented experimentally (Harrison and Boyett, 1995), and was found necessary to reproduce the experimentally observed rise of diastolic calcium at high pacing rates. The second is the compartmentation of the myoplasm into a submembrane space near the sarcolemma and the rest of the myoplasm. This distinction recognizes the fact that the concentration in the submembrane space, which influences both the NaCa exchange current and the calcium-dependent inactivation of ICa, is much larger during release than the average bulk concentration inside the myoplasm, due to the proximity of this space to dyadic junctions. Finally, consistent with the discrete picture of sparking events, the model builds in a diffusional delay between the network SR and the JSR that distinguishes between the average JSR concentration of recruited and unrecruited dyadic junctions. We shall show here that this distinction leads to the interesting effect that this delay can enhance SR release, as opposed to always reducing it, as in previous models where the JSR is treated as a single pool (Jafri et al., 1998; Glukhovsky et al., 1998; Luo and Rudy, 1994).
We explore the dynamics of calcium cycling by pacing the model using a clamped AP waveform that is fitted to those used in the experiments of Chudin et al. (1999). At high pacing rates, the model exhibits sustained calcium alternans, consistent with the experimental findings. Furthermore, we vary parameters of the model to identify which parts of the complex calcium cycling machinery are essential for alternans and which ones have a secondary role. We identify two essential instability mechanisms for alternans. The first is a steep relation between release and SR load, and the second is a steep relation between release and diastolic calcium that is due to calcium-induced inactivation of ICa. In addition, we find that the calcium dynamics can become chaotic even when paced with a periodic AP clamp.
To gain further insight into both alternans and chaotic dynamics, we reduce the ordinary differential equations (ODEs) of the model to a two-dimensional discrete map that relates the SR and cytosolic calcium concentrations at one beat to those at the previous beat. This map allows us to obtain a stability condition for the onset of calcium alternans with AP clamp pacing in terms of the slopes of the relations between release and SR load, and release and diastolic calcium, as well as the strength of the calcium uptake into the SR. In addition, this map allows us to show that the existence of chaos is due to a nonmonotonic relation between the peak calcium transient at one beat to that at the previous beat.
This article is organized as follows. In the next section, we describe the various intracellular compartments used in the model, and define the corresponding calcium concentrations. In the section "SR Release Flux", we describe a model for the SR release current. The calcium dependent membrane currents are then described in the section, "Calcium Fluxes across the Sarcolemma". In the section "Dynamics of Calcium Cycling", we integrate these currents into a set ODEs that are used to study calcium cycling. In the section "Voltage Clamp Pacing", we present the results of the model when paced under an AP clamp protocol. The reduction of this model to a map is then described in the section "Analysis of Beat-to-Beat Dynamics", and the results are used to analyze the dynamical behavior of the model. The results are discussed in the "Discussion" section, followed by "Conclusions".
|
INTRACELLULAR COMPARTMENTS AND CONCENTRATIONS
|
|---|
A schematic diagram of the intracellular compartments relevant to calcium cycling is given in Fig. 1. The SR is a spatially diffuse network of tubules and cisternae that is composed of two distinct parts. The first is the network SR (NSR), which is a meshwork of tubules that enwraps the myofilaments. These tubules branch out from the NSR into flattened elliptical sacs, referred to as the junctional SR, that position themselves close to the surface of the cell sarcolemma (Forbes et al., 1985; Franzini-Armstrong et al., 1999). The volume in between the JSR and the sarcolemma is the dyadic junction, wherein RyR channels in the JSR membrane are in close proximity to L-type channels in the sarcolemma. In ventricular myocytes, the sarcolemma forms a uniform array of deep invaginations into the cell, referred to as T-tubules (Soeller and Cannell, 1999), which effectively distribute dyadic junctions uniformly throughout the cell.
The flux of calcium from the SR to the myoplasm, shortly after a voltage depolarization, is due to a summation of several thousand discrete local release events (Cannell et al., 1994; Cleemann et al., 1998). During the release process, calcium in the SR is being drained into the myoplasm via JSR compartments that are being depleted by local calcium release. To describe the calcium concentration within the SR during this release process, we first assume that the SR is composed of Njsr identical JSR compartments, each with a volume, 
, along with the bulk NSR with volume 
. This gives a total SR volume of 
. The total calcium concentration within the kth JSR compartment will be denoted by 
and the concentration in the NSR is cnsr. At any given time, the JSR compartments can be divided into those that are being drained due to a spark in the local dyadic space, and those that are not. The average total calcium concentration within JSR compartments that are not being drained at time t will be denoted by 
. This average is simply given by
 | (1) |
where the summation is over the set of NU(t) unrecruited JSR compartments. The average total calcium concentration within the entire SR network, which includes all the JSR and the NSR, will be denoted by cj(t) and is given by
 | (2) |
where the summation is over both recruited and unrecruited JSR volumes. It is important to distinguish between cj and 
, since in general cj 
. This is due to the finite time necessary for to relax to cj via the diffusion of calcium ions within the SR network. During a local release, at the kth dyadic junction, for instance, 
will drop below cj as the JSR compartment is depleted. However, can be larger than cj since the set of unrecruited JSR compartments are not depleted. This feature is quite different from previous models (Jafri et al., 1998; Glukhovsky et al., 1998; Luo and Rudy, 1994), where the entire JSR is treated as a single pool that is being replenished by the NSR. In such models, the JSR concentration is always less than the SR concentration during release.
To describe the calcium outside of the SR, we first divide this space into the bulk myoplasm (volume 
i, average free calcium concentration ci), and a submembrane space (volume 
, average free calcium concentration cs), which corresponds to the space in between the bulk myoplasm and the sarcolemma (see Fig. 1). To clarify further, cs is the average concentration that calcium-dependent membrane channels sense. Since membrane currents are delivered into , and since 
(in this article we use 
10), the concentration changes in the submembrane space are much larger and faster than those in the bulk myoplasm. This effect can potentially influence calcium-dependent ion channels at the sarcolemma. Also, the free calcium concentration within the kth dyadic junction will be denoted by 
It is important to note that the average volume of a dyadic space is 10-9 times smaller than that of the cell volume. Thus, during a calcium spark in the kth dyadic junction the concentration 
can rise rapidly to levels >100 µM (Peskoff and Langer, 1998), whereas the concentration in the submembrane space will rise to 5–10 µM.
|
SR RELEASE FLUX
|
|---|
Spark recruitment rate
To model the release process at the whole cell level, it is first necessary to determine the relationship between the number of release events (sparks) and the corresponding voltage stimulus. Since RyR channels are triggered by nearby L-type channels, it is reasonable to expect that spark recruitment is related to the total calcium current entering the cell via the L-type channels (ICa). Collier et al. (1999) have measured whole cell ICa simultaneously with calcium spark occurrence in rat ventricular myocytes, albeit under experimental conditions of reduced ICa, which allowed for the visualization of individual sparks. They showed that the time course (i.e., time constant of decreasing occurrence) of calcium sparks with sustained depolarization was statistically identical to the time course (i.e., the time constant of inactivation) of the ICa current. This confirmed the conclusions of an earlier study that was similar except that ICa was not measured (Lopez-Lopez, 1995). A more recent article (Sah et al., 2002) again shows in a quantitative manner that the time course of spark occurrence parallels the time course of ICa, using novel voltage clamp protocols, such as triangular and action potential waveforms.
To quantify the relationship between spark occurrence and the ICa current, we will denote the number of sparks recruited at time t to be N(t). Then the rate of spark recruitment is dN(t)/dt, which is the number of sparks recruited in a whole cell per unit time. This rate can be computed directly from experimental measurements of spark occurrence. Using the experimental data presented in the work of Collier et al. (1999), we compared the rate of spark occurrence and the ICa current for a depolarization to a holding potential of 30 mV (see Fig. 4 in their article). In Fig. 2, we plot the computed spark rate as a function of whole cell ICa. It is clear from the plot that, under the physiological conditions of the measurement and for the given depolarization potential, the measured spark rate is to a good approximation proportional to whole cell ICa, i.e., dN(t)/dt 
ICa(t).
View larger version (8K):
[in this window]
[in a new window]
|
FIGURE 4 (A) Plot of the ratio x = APD/T versus T. The circles correspond to values of x computed from the experimental AP clamps, and the line is the corresponding fit. (B) Maximum and minimum of bulk myoplasmic calcium concentration during steady state as a function of pacing period T at relatively slow pacing rates. The solid lines are for the case when the internal sodium concentration increases with decreasing period according to Eq. 17. The filled circles correspond to the experimental data points from Chudin et al. (1999). The dashed line corresponds to the case when intracellular sodium is fixed at Nai = 10 mM.
|
|
View larger version (10K):
[in this window]
[in a new window]
|
FIGURE 2 Spark rate versus whole cell ICa using the experimental data of Collier et al. (1999). The straight line corresponds to a linear fit.
|
|
It is known experimentally that when the calcium content of the SR is increased, the frequency of spontaneous sparks in a resting myocyte also increases (Cheng et al., 1993; Lukyanenko et al., 1996, 2000). This dependence between spark occurrence and SR content implies that RyR channels are sensitive to the calcium concentration within the local JSR compartment. Now, since JSR compartments, which already have been depleted due to a spark, probably cannot be recruited until they have had enough time to refill, we expect that the rate of spark recruitment should depend on the average calcium concentration within unrecruited JSR compartments (). Thus, we model the JSR calcium dependence of the whole cell spark rate using
 | (3) |
where the function A() gives the JSR load dependence, and where g is a proportionality constant.
Calcium release during a spark
The local release flux during a spark will be dictated by the gating kinetics of the RyR cluster and the calcium gradients in the dyadic space. However, because the detailed properties of a cluster of RyR channels are not well known, we will describe local release using a simple phenomenological model based on very general considerations. First, we will assume that a spark that is activated at a given time t', in the kth dyadic junction, is determined primarily by the local JSR concentration at the time of activation 
rather than the concentration in the local dyadic space 
This is a reasonable assumption because the free calcium concentration gradient across the RyR channels is extremely large, and once CICR is induced, the subsequent flux will rise to a value that is determined by the amount of releasable calcium in the JSR. We will also assume that once a spark is recruited it will have a constant lifetime denoted by 
. This assumption is motivated by fluorescence measurements of calcium sparks (Cleemann et al., 1998; Niggli, 1999), which show that sparks have a well defined lifetime around 10–30 ms. Given these simplifying assumptions, the flux at time t, from the kth JSR compartment, during a spark that is activated at time t' 
t, can be written as
 | (4) |
where the function Jk gives the JSR load dependence of the local flux amplitude at the kth dyadic junction.
Once the whole cell ICa is activated, sparks are stochastically recruited at many dyadic junctions distributed throughout the cell. Thus, at a given time t', an ensemble of release fluxes given by Eq. 4 is recruited. Here, we will argue that this ensemble of spark fluxes has well-defined averaged properties. This averaging procedure will allow us to drop the superscript k in Eq. 4, and will simplify the subsequent analysis of spark summation. First, we note that 
in Eq. 4 denotes the JSR concentration at the time t' when a spark is just being recruited. Thus, the kth JSR compartment has not yet been depleted and will have a concentration that is roughly the same as the average concentration in unrecruited JSR compartments (t'). Hence, we can make the approximation that 
Second, we assume that the lifetime and amplitude of a spark in the kth dyadic junction is roughly the same as the average over the ensemble of recruited sparks. To be more precise, if we define the average 
and 
where the brackets denote an average over the ensemble of recruited sparks, then 
and Jk 
J. This claim is equivalent to the statement that the distribution of spark lifetimes and amplitudes is narrowly peaked around their averages. This assumption is motivated by fluorescence images of calcium sparks that show that sparks in different regions of the cell essentially have the same spatial and temporal properties (Cleemann et al., 1998; Niggli, 1999). Thus, the local flux during a spark in the kth junction 
can be represented by an average spark flux 
where 
and 
To asses the validity of Eq. 4, we have also developed a set of ODEs that describe the calcium dynamics within a dyadic junction. This model, given in detail in Appendix A, simulates a calcium spark using a simple CICR gating kinetics that is initiated by a short pulse of calcium injected into a simulated dyadic space. Local flux termination is then modeled using a combination of diffusion-driven deactivation and local JSR depletion. In Appendix A, we show that Eq. 4 approximates well the local current due to a simulated local release. However, it is important to note that Eq. 4 does not rely on any specific mechanism underlying the local calcium dynamics, and should hold under more general conditions. Thus, our simple model of the calcium dynamics at the dyadic junction level serves primarily to validate Eq. 4 in a particular case, but is not central to the overall implementation of the model.
Summation of calcium sparks
The total flux draining the SR at time t, which will be denoted by Ir(t), is given by the summation of local fluxes due to sparks recruited at times t' t. To compute this sum, we divide the interval of time 0 t' t into M bins of duration 
t', such that t = Mt'. Consequently, the flux draining the SR at time t = Mt' is given by the discrete sum
 | (5) |
where N(it') is the number of sparks recruited during the time interval (i - 1)t' t' it'. If we now take the continuum limit 
and use the expression for the average spark flux Ispark(t, t') (Eq. 4 where the superscript k is dropped), the above discrete sum becomes an integral
 | (6) |
where t' = 0 is a time origin when Ir(0) = 0. Finally, differentiating both sides of this equation with respect to t, we obtain
 | (7) |
which will be our basic equation for the release flux Ir(t). Using Eq. 3, the first term on the right-hand side of Eq. 7 can be written as 
which reveals that the JSR load dependence of release is governed by the function 
SR load dependence of release
The release of calcium as a function of SR load has been addressed experimentally (Bassani et al., 1995; Shannon et al., 2000b). In the experiment by Bassani et al. (1995), the fraction of SR calcium released as a function of SR load at fixed trigger ICa was measured in ferret ventricular myocytes. The authors found that for low SR loads [Ca]SRT < 50 µmol/l cytosol, almost no calcium is released. At normal loading [Ca]SRT 90 µmol/l cytosol, about half of the SR is released, whereas for slightly larger loading [Ca]SRT 100 µmol/l cytosol, almost 70% of the SR is released. Here, [Ca]SRT denotes the total calcium concentration in the SR, which is the quantity measured in the experiment. An important finding of this experiment is that the fraction of calcium released from the SR depends on the initial calcium concentration of the SR in a highly nonlinear fashion. To model this dependence of release on SR load, we will pick a functional form for Q(), which will reproduce the essential features of the experimental results. Since the dependence of release on load only enters through the function Q(), it is clear that this function itself should depend on in a nonlinear fashion. There can be several explanations for such a nonlinear dependence of local release on JSR load. One possible scenario is that by virtue of the gating kinetics of the local RyR channels, the amount of calcium released at a single dyad J(), may change nonlinearly with JSR load. Another independent possibility is that the sensitivity of the RyR channels may depend on in a nonlinear fashion. In this case, it would be the load dependence of the spark rate A() that is a nonlinear function. However, the advantage of our phenomenological approach is that the function Q() can be chosen by appealing directly to experimental data, without a complete model of the SR load dependence of release.
JSR-SR diffusional relaxation
To complete the model, we take into account the fact that the average total calcium concentration in the SR (cj) is not the same as the average concentration of unrecruited JSR compartments (). This is due to the finite time necessary for calcium to equilibrate, via diffusion, over the entire volume of the SR. This effect can be incorporated into the model by letting relax to cj using 
= (cj - )/
a, where a is a relaxation time. A rough estimate of a, which has not yet been measured experimentally, is the average time for a calcium ion to diffuse from a local JSR compartment to the bulk NSR. If we take the effective diffusion coefficient of calcium ions in the SR to be 150 (µm)2/s (Stern et al., 1999), then the time for an unobstructed calcium ion to diffuse over a distance 1 µm, (which is roughly the average distance between the NSR and the sarcolemma) is 7 ms. However, the calcium ions are confined within the SR, which is a complex tubular network of diameter 25–50 nm (Forbes et al., 1985), and so the diffusion (translocation) time can actually be much longer. Thus, we will consider relaxation times within a fairly wide range of 1–100 ms, and discuss in detail the role of this time constant on the dynamical behavior of the calcium system.
|
CALCIUM FLUXES ACROSS THE SARCOLEMMA
|
|---|
L-type calcium current
In this model, we have approximated the rate of spark recruitment to be proportional to the whole cell ICa current entering the cell at that time. Hence, a correct formulation of whole cell ICa is an essential component of the model. The whole cell ICa current can be written as ICa = MPoiCa, where M is the number of L-type channels in the cell, Po is the time-dependent open probability of a single channel, and iCa is the single-channel current. We model the open probability of the L-type channel using Po = d
fq, where d is an instantaneous voltage dependent activation gate, f is a slow voltage dependent inactivation gate variable, and where q describes calcium-induced inactivation. It is important to note that L-type channels, within dyadic junctions where a spark has just been activated, will see a calcium concentration that is much larger than cs. Thus, these channels will have a much lower open probability as compared to the open probability of channels in unrecruited dyadic junctions. Hence, at a given time t, calcium entry into the whole cell via ICa(t) should be dominated by those L-type channels within unrecruited dyadic junctions. Since the average concentration sensed by these channels is roughly the same as cs, then the calcium-dependent gate variable q should depend on cs.
The kinetics of calcium inactivation is modeled using a simple first order scheme that yields
 | (8) |
The detailed parameters of the gate variables are given in Appendix B, and are chosen so that the time course of ICa is qualitatively similar to the most recent experimental measurements in rabbit myocytes (Puglisi et al., 1999).
To explore the role of calcium-induced inactivation on the overall dynamics of the system, it is essential to also investigate more general inactivation schemes. However, we found that simply generalizing Eq. 8 to higher order kinetics led to an ICa time course that was different from experimentally measured currents (Linz and Meyer, 1998; Puglisi et al., 1999). Following an approach similar to Fox et al. (2002), we modeled calcium-induced inactivation using a simple phenomenological scheme with
 | (9) |
and with q a constant in the range 10–50 ms. With this model formulation the exponent 
can be varied, and a steeper cs dependence of calcium-induced inactivation can be explored.
NaCa exchange current
The NaCa exchange current is taken directly from Luo-Rudy II (Eq. 31 in Appendix B). The exchange current is sensitive to the intracellular sodium concentration (Nai). Higher Nai reduces the transmembrane sodium gradient and so reduces the efficiency of the NaCa exchange current at extruding calcium from the cell, which should lead to higher intracellular calcium. An investigation by Harrison and Boyett (1995) used guinea pig ventricular myocytes to study pacing rate dependent increases in twitch shortening. They showed that when the pacing rate was increased from 0.5 to 3 Hz, twitch shortening increased by 34%, whereas if sodium entry (INa) was blocked by tetrodotoxin, twitch shortening increased by only 8%. This result demonstrates that rate dependent changes in Nai could have a substantial effect on intracellular calcium. Thus, we felt it imperative to include a phenomenological description of Nai, where Nai is taken to be an increasing function of the pacing rate.
Role of the submembrane space
As mentioned in "Intracellular Compartments and Concentrations", we have divided the myoplasm into two spaces: a submembrane space, which is the space in the vicinity of the sarcolemma, and the bulk myoplasm, which is the remainder of the cell volume. This compartmentation of the myoplasm is motivated by the recent experimental and modeling work of Weber et al. (2001). Here, it was shown in particular that the NaCa exchange current INaCa, due to its proximity to dyadic junctions, senses a calcium concentration that is different from the global cytosolic calcium concentration. Consequently, modeling studies showed that the time course of INaCa served as a much more efficient calcium efflux mechanism when it is coupled to the larger and faster concentration changes at the submembrane space (cs) than when it was coupled to the bulk concentration (ci). Likewise, by virtue of calcium-induced inactivation, the time course of ICa should depend on the calcium concentration in the vicinity of the sarcolemma. This effect is particularly important since the calcium released from RyR channels during calcium release is delivered into dyadic junctions that are close to the sarcolemma. Thus, L-type channels sense a much more rapid rise in calcium than that given by the rise of the global cytosolic calcium concentration. To model the effect of a submembrane space, we let the membrane currents ICa and INaCa flow into a compartment of volume 
before diffusing to the bulk myoplasm via a simple relaxational current Id = (cs - ci)/s, with s a relaxation time that we will take to be in the range 5–10 ms.
|
DYNAMICS OF CALCIUM CYCLING
|
|---|
In this section, we study calcium cycling by incorporating the various currents described in previous sections into a single model. Since our formulation of the SR release flux does not explicitly depend on the calcium concentration gradient between the JSR and the submembrane space, we are free to formulate the dynamics of the SR calcium concentration in terms of the total concentration in the SR, rather than the free calcium concentration. An advantage of using the total SR concentration is that it is the quantity that is measured directly in experiments (Shannon et al., 2000b). Calcium buffering in the cytosol will be taken into account by incorporating the buffering to troponin C, SR, and calmodulin sites, since these are the sites that bind the majority of calcium released from the SR (Shannon et al., 2000a). The buffering to SR and calmodulin sites are fast and will be treated as instantaneous, whereas the time dependent kinetics for the buffering to troponin C will be accounted for.
The equations for calcium cycling can be written as
 | (10) |
where the flux out of the SR satisfies
 | (11) |
and where
 | (12) |
Here, 
and 
describe the time dependent buffering to troponin C in the submembrane and bulk myoplasm respectively. The functions ß(cs) and ß(ci) account for the instantaneous buffering to SR and calmodulin sites. The details of the buffering kinetics and parameters are given in Appendix C. Also, V denotes a time-dependent AP waveform.
The free concentration in the cytosol (ci) and the submembrane space (cs) are in units of µM, whereas the total concentration in the SR and the JSR is scaled by a constant factor 
, so that cj is in units of µmol/l cytosol. All fluxes have been divided by the volume of the cytosol 
, and have units of µM/s. Thus, currents can be converted to amperes by multiplying by the factor 2
, where F is Faraday's constant. Also, the volume of the cytosol will be taken to be 10 times the volume of the submembrane space (
). The parameters used in the model are given in Tables 1–4. It should be noted that the above equations do not mathematically forbid unphysiological draining of the SR, where cj becomes less than zero. However, for physiologically plausible parameters, such as the ones used in our simulations, cj is always positive.
The constant g will be taken to have units of sparks/µM, so that gICa gives the rate, in units of sparks/s, that sparks are recruited in the whole cell. The release function Q() is taken to have the simple form
 | (13) |
where u is an adjustable constant, and where s is fixed by the condition that the release function is continuous. Here, the flux Q() is expressed in units of 10-3 µM/s, which, for a typical JSR load of = 100 µM, yields a local release current of 1 pA, which is roughly the same as the experimental estimate of local release flux (Cheng et al., 1993). With the above choice of units we will use g = 1.5 x 104 sparks/µM. Thus, for a total calcium entry, during a 1-s cycle, of 
4 x 104 sparks are recruited.
The constant u is adjusted so that the fractional release is consistent with the experimental results of Bassani et al. (1995). The fractional release is defined as
 | (14) |
where 
is the total SR concentration just before an AP upstroke, and where 
is the minimum total SR concentration shortly after the upstroke. To compute the fractional release, we clamp the voltage from Vmin = -80 mV to Vmax = 30 mV. Under these conditions, ICa quickly peaks and decays during which the SR concentration crosses a minimum value. By computing the minimum value for various initial SR loads, we can compute the function f(). In Fig. 3 A we plot the fractional release f versus the initial SR load for the parameters given in Tables 1–4.
An important experimental observation that has to be incorporated into any model of EC coupling is graded release. In an important experiment, Wier et al. (1994) measured peak SR flux and peak ICa during depolarizations to various test potentials. In a similar manner, we held the membrane voltage at a V(t) of -80 mV, and depolarized to various test potentials (Vmax). As expected, and similar to Wier's results, peak ICa and peak Ir were bell shaped as a function of Vmax (Fig. 3 B), and peak SR flux was graded with respect to peak ICa. The parameters were adjusted to give a gain 
of 10 for Vmax in the range 0–30 mV (Wier et al., 1994).
|
Voltage Clamp Pacing
|
|---|
Pacing protocol
In the experiment of Chudin et al. (1999), intracellular calcium transients were measured as the cell was paced with various AP clamps. The AP clamp used in that experiment can be modeled effectively using
 | (15) |
where T is the period (pacing cycle length), x = APD/T, where APD denotes the action potential duration, and where m denotes the mth beat. The parameter x fixes the shape of the AP, and depends on the pacing period. The action potentials are fitted using Vmin = -80 mV and Vmax = 30 mV. The ratio x = APD/T was fitted to the experimental values using a function of the form
 | (16) |
In Fig. 4 A, we plot the experimental values for x as a function of T, as well as the values given by a fit with a = 2/3.
Intracellular sodium and the rise of the calcium transients
The average whole cell sodium concentration Nai depends on the amount of sodium injected into the cell during pacing, and so depends on the pacing period T. Since our model does not describe intracellular sodium dynamics, we introduce a function Nai(T), which gives the internal sodium concentration at steady state as a function of period T. Faber and Rudy (2000) have suggested that Nai increases from 10 mM to close to 15 mM, as the pacing period is decreased from 2 s to 100 ms. Hence, we model the intracellular sodium concentration using a simple function of the form
 | (17) |
where the constants a, b are chosen so that Nai(1 s) = 9.6 mM and Nai(.2 s) = 18.5 mM. In Fig. 4 B, we plot 
and 
after steady state has been reached for the above set of parameters. On the same graph we include the experimental data of Chudin et al. (1999). As we can see, both quantities rise with decreasing period, as - remains fairly constant. On the same graph we plot the same quantities when intracellular sodium concentration is kept constant at Nai = 10 mM. In this case there is only a small rise in . Hence, in this model, the overall rise in calcium concentration during rapid pacing is primarily due to the accumulation of sodium in the cell.
Calcium transients and currents during rapid pacing
We studied the effects of pacing rate by pacing the calcium system with the clamped AP waveform at different cycle lengths (T). Results of pacing the model at a slow rate, with pacing cycle length of T = 1 s, and a fast rate, with cycle length T = 0.25 s, are shown in Figs. 5 and 6, respectively. Both sets of figures show the steady state values of ci, whole cell ICa, and the NaCa exchange current INaCa. The results of Fig. 6 show clearly that as the whole cell calcium transient alternates, so do the calcium-dependent currents. We explored the dynamics of the model more fully by calculating the peak calcium transient values at different pacing periods (Fig. 7). On the same graph, we plot the experimental data of Chudin et al. (1999). It is clear from the graph that as the period is decreased, the calcium system undergoes a period-doubling bifurcation, qualitatively similar to the experimental findings.
View larger version (16K):
[in this window]
[in a new window]
|
FIGURE 5 Calcium transients and currents during steady state for a pacing period of T = 1 s. (A) Plot of ci versus t. (B) Plot of ICa versus t. (C) Plot of INaCa versus t.
|
|
View larger version (21K):
[in this window]
[in a new window]
|
FIGURE 6 Calcium transients and currents during steady state for a pacing period of T = 0.250 s. (A) Plot of ci versus t. (B) Plot of ICa versus t. (C) Plot of INaCa versus t.
|
|
View larger version (12K):
[in this window]
[in a new window]
|
FIGURE 7 Peak values of the bulk myoplasmic calcium concentration during steady state as a function of pacing period T at fast pacing rates. The open circles correspond to the experimental results of Chudin et al. (1999). The filled circles correspond to the model predictions.
|
|
|
ANALYSIS OF BEAT-TO-BEAT DYNAMICS
|
|---|
Map reduction with linear instantaneous buffering
To understand the nonlinear dynamics of the calcium system, it is useful to reduce the dynamics of the model ODEs to a system of discrete maps. Using this approach, it is possible to understand the dynamical instability that is responsible for calcium alternans. To make the analysis more tractable, we first focus on the simple case where all cytosolic calcium buffering is linear and instantaneous; in this case, ß(cs) = ßs, and ß(ci) = ßi, where ßi and ßs are constants. This simplification is introduced mainly because it allows us to derive a simple analytical condition for the onset of calcium alternans, which makes transparent the physiological mechanisms that underlie this phenomenon. In the section "Map reduction with instantaneous nonlinear buffering", we then show that a more elaborate map derived with instantaneous, albeit nonlinear, buffering can yield predictions of calcium dynamics that are in reasonable quantitative agreement with the full model ODEs.
Recall that our basic model (Eqs. 10–12) consists of four variables that represent the average concentrations in the cytosol, ci(t); the submembrane space, cs(t); the JSR, (t); and the SR, cj(t). Furthermore, let us define 
ci(nT), 
cs(nT), 
cj(nT), and 
(nT), to be these concentrations at a time t = nT just before the nth AP upstroke. The abruptness of the AP upstroke makes this a natural choice of variables for the discrete map. A mapping between these concentrations at the nth and nth + 1 AP upstrokes can be constructed by integrating the model ODEs from time t = nT to time t = (n + 1)T. A major simplification is that the concentrations in the cytosol and the submembrane space are to a very good approximation equal before an AP upstroke (
) and only differ significantly from each other during the upstroke of the calcium transient. This is due to the fact that the volume of the submembrane space is much smaller than the volume of the cytosol 
and hence cs relaxes quickly to ci. Hence, the two concentrations are essentially equal except when ci is rapidly varying. Similarly, the concentrations in the SR and in the JSR can be assumed to be equal just before an AP upstroke (
) as long as the time for calcium to diffuse between these two compartments (a 50 ms) is much smaller than the pacing interval 
Using the fact that both and , it is straightforward to obtain from Eq. 10 the map
 | (18) |
 | (19) |
where 
= ßi/(1 + ßi
s/ßsi), and where we have defined the quantities
 | (20) |
which represent the total amount of calcium release (R) and uptake (U) in one beat as well as the net total amount of calcium () that enters the cell. Both R and U are always positive whereas can be either positive or negative depending on whether the amount of calcium that enters the cell through the L-type channels is larger or smaller than the amount of calcium extruded from the cell by the NaCa exchange current in one beat. The currents that appear in the integrals above must be calculated by integrating the full set of ODEs (Eqs. 10–12) for the variables (cs(t), ci(t), cj(t), and (t)) using as initial conditions cs(nT) = ci(nT) = and cj(nT) = (nT) = , and values of the gate variables d(nT) and f(nT) for ICa that are determined by the periodic clamped AP waveform.
Calcium alternans
In the rest of this section, we further analyze this map to pinpoint the condition for the onset of alternans. The fixed point of the map = 
and = is defined by the conditions R = U and = 0, which implies that at steady state the total uptake during one cycle should be equal to the total release, and the net total flux into the cell over one period should be zero.
To obtain a simple analytical condition for the onset of alternans, we assume that the beat-to-beat variation of the net free calcium that enters the cell in one beat is much smaller than the variation of the total amount of calcium pumped into the SR (
). This limit is physiologically relevant since the total calcium uptake in one beat is typically 3–4 times larger than the amount of calcium extruded from the cell via the NaCa exchange current (except in diseased states such as heart failure where the uptake and extruded amounts can become comparable). In this limit, the total diastolic amount of calcium in the cell (i.e., including both the myoplasm and the SR) can be assumed to remain approximately constant from beat to beat, or
 | (21) |
The above relation allows us to reduce the system of two coupled m
TOP
ABSTRACT
INTRODUCTION
(Vmax), and the dashed line is 
(Vmax). Both curves have been normalized to their peak values. The initial conditions are identical to those given above except that the initial total SR load is fixed at cj(0) = 100 µmol/l cytosol.