- Open Access
Modeling CICR in rat ventricular myocytes: voltage clamp studies
Theoretical Biology and Medical Modelling volume 7, Article number: 43 (2010)
The past thirty-five years have seen an intense search for the molecular mechanisms underlying calcium-induced calcium-release (CICR) in cardiac myocytes, with voltage clamp (VC) studies being the leading tool employed. Several VC protocols including lowering of extracellular calcium to affect Ca2+ loading of the sarcoplasmic reticulum (SR), and administration of blockers caffeine and thapsigargin have been utilized to probe the phenomena surrounding SR Ca2+ release. Here, we develop a deterministic mathematical model of a rat ventricular myocyte under VC conditions, to better understand mechanisms underlying the response of an isolated cell to calcium perturbation. Motivation for the study was to pinpoint key control variables influencing CICR and examine the role of CICR in the context of a physiological control system regulating cytosolic Ca2+ concentration ([Ca2+] myo ).
The cell model consists of an electrical-equivalent model for the cell membrane and a fluid-compartment model describing the flux of ionic species between the extracellular and several intracellular compartments (cell cytosol, SR and the dyadic coupling unit (DCU), in which resides the mechanistic basis of CICR). The DCU is described as a controller-actuator mechanism, internally stabilized by negative feedback control of the unit's two diametrically-opposed Ca2+ channels (trigger-channel and release-channel). It releases Ca2+ flux into the cyto-plasm and is in turn enclosed within a negative feedback loop involving the SERCA pump, regulating[Ca2+] myo .
Our model reproduces measured VC data published by several laboratories, and generates graded Ca2+ release at high Ca2+ gain in a homeostatically-controlled environment where [Ca2+] myo is precisely regulated. We elucidate the importance of the DCU elements in this process, particularly the role of the ryanodine receptor in controlling SR Ca2+ release, its activation by trigger Ca2+, and its refractory characteristics mediated by the luminal SR Ca2+ sensor. Proper functioning of the DCU, sodium-calcium exchangers and SERCA pump are important in achieving negative feedback control and hence Ca2+ homeostasis.
We examine the role of the above Ca2+ regulating mechanisms in handling various types of induced disturbances in Ca2+ levels by quantifying cellular Ca2+ balance. Our model provides biophysically-based explanations of phenomena associated with CICR generating useful and testable hypotheses.
Contraction of cardiac muscle is triggered by a transient rise in intracellular Ca2+ concentration [Ca2+] myo . Sarcolemmal (SL) membrane depolarization triggers Ca2+ influx from the extracellular medium by opening dihydropyridine (DHP)-sensitive L-type Ca2+ channels. Following diffusion across a small sub-membrane dyadic space, this influx activates ryanodine receptors (RyRs) controlling ryanodine-sensitive Ca2+ release channels in the junctional portion of the sarcoplasmic reticulum (jSR). Fabiato and Fabiato  named the process calcium-induced calcium release (CICR). Ca2+ subsequently diffuses from the dyadic space into the cytosol. Ultimately, intracellular Ca2+ concentration [Ca2+] myo is returned to resting levels by combination of: (a) Ca2+ buffering in the dyadic space and cytosol; (b) sequestration of Ca2+ by sarcoplasmic/endoplasmic reticulum Ca2+-ATPase (SERCA)-type calcium pumps lining the longitudinal portion of the sarcoplasmic reticulum (LSR); and (c) Ca2+ extrusion from the cytosol by Na+/Ca2+ exchangers and Ca2+-ATPase pumps on the sarcolemmal membrane.
CICR in cardiac muscle exhibits both graded behavior and a high gain. Graded behavior refers to the obser-vation that SR Ca2+ release is proportional to the influx of trigger Ca2+ , whereas high gain indicates that the SL trigger current elicits a high SR Ca2+ release flux. Graded Ca2+ release with high gain is somewhat paradoxical according to Stern , in that the positive feedback inherent in such high-gain systems tend to produce regenerative, nearly all-or-none release rather than graded release. Several deterministic models have been developed to explain excitation-contraction (E-C) coupling [4, 5], but none of them can explain the mechanism of graded release at high gain over a wide range of values for sarcolemmal Ca2+ current. Stern  proposed that such a gradation paradox might be explained if the stimulus for Ca2+ release by RyRs were actually the local nanodomains of [Ca2+] generated by nearby L-type channels, rather than the global cytosolic [Ca2+] myo . According to this hypothesis, graded control of macroscopic SR Ca2+ release can be achieved by graded statistical recruitment of individual, autonomous, all-or-none stochastic release events . In these studies, a distributed differential model of high order that included dynamic interactions between large numbers of individual channels was used to demonstrate this concept. However, rather large amounts of computation time are required with distributed stochastic models of this type. Additional models have sought to characterize the Ca2+ release complex, including several [7–9] based on the stochastic release process adopted by Stern et al. These statistical models have solved the graded release problem, however, they too are complicated and computationally very expensive. Other models based on the simplified local control model of CICR developed by Hinch et al.  sought to adopt a lower order description of the E-C coupling process [11, 12] by making an approximation of rapid equilibrium in the dyadic space. The latency from onset of Ca2+ entry via the I Ca,L channel to triggered SR Ca2+ release is known to increase with decrease in the magnitude of I Ca,L,TT , the modeling of which is made possible by considering Ca2+ diffusion in the dyadic medium. These models [11, 12] also approximate the SR as a single volume compartment with no distinction between junctional versus the longitudinal (network) SR compartments. However, recent work  points towards the important role of the Ca2+ refilling rate from the network to junctional SR in controlling RyR release termination via the luminal sensor. Shiferaw et al.  developed a computationally tractable model of Ca2+ cycling to represent the release of calcium from the SR as a sum of spatially localized events that correspond to Ca2+ sparks, assuming the recruitment rate of Ca2+ sparks is directly proportional to the whole-cell I Ca,L current. This assumption overlooks the complex calmodulin mediated interaction (calcium dependent facilitation (CDF) and calcium dependent inactivation (CDI)) of the I Ca,L channel with calcium in its vicinity. It also demands a large amount of computation.
Numerically, distributed as well as statistical models tend to be computationally expensive due to the in-herent repetition involved in the computation. In a spatially distributed model, simultaneous solution for dynamics in identical compartments distributed in space would amount to a large computational cost. In statistical models inference is drawn based on multiple runs of identical events which translate into a prolonged simulation time. Hence, these models are cumbersome to implement, particularly in larger multiple-cell simulations. Consequently, we consider a deterministic approach to the characterization of CICR. Specifically, we develop a lumped model of the Ca2+ release complex that includes: (a) a sub-sarcolemmal dyadic cleft space separating the SL and jSR membranes; (b) a single DHP-sensitive Ca2+ channel on the SL membrane; and (c) a single equivalent Ry-sensitive channel arranged symmetrically on the opposing jSR membrane that represents the output of a local cluster of Ry-sensitive channels facing the DHP-sensitive channel. Based on morphological data compiled by Bers , we further assume that each ventricular cell contains 10,000 of these dyadic Ca2+ release units, and that they are associated with the fraction of the SL membrane that is coupled with the jSR. That is, we partition the sarcolemma into free and dyadically coupled SL membrane, and associate each with a different fluid compartment: the cell cytosolic medium in the case of the free SL membrane, and the dyadic cleft space medium in the case of the dyadic-coupled fraction. In a sense, we build on Stern's  local domain concept by considering the aforementioned local nanodomains identical, but focusing on the nonlinear dynamics of the two different types of Ca2+ channel in the dyadic coupling unit. Our deterministic model although is very descriptive, is computationally tractable and has a run time of 21 sec (including recording of 73 variables of type double on a data file) for 1 cycle of 4 Hz voltage clamp stimulation.
Rat ventricular myocytes were prepared from 200-300 g male Sprague Dawley rats by dissociation with collagenase, as previously described . All experiments were performed under conventional whole cell recording conditions with a List EPC-7 patch clamp, recording fluorescence from nearly the entire cell, as described by Fan and Palade . Recordings from an individual cell were rarely extended beyond 10 min in order to reduce as much as possible both escape of dye from the cell and Ca2+ current rundown. External solution in the bath was normal Tyrode (1 mM Ca2+) with Cs+ substituted for K+ for purposes of blocking inward rectifier K+ currents. The internal solution in the pipette contained Cs aspartate supplemented with 20 mM CsCl, 3 mM Na 2 ATP , 3.5 mM MgCl 2 and 5 mM HEPES. Holding potential used was -40 mV.
All simulations and analysis were performed on a 2.8 GHz Intel® Core™2 Duo CPU-based computer using Microsoft Windows XP operating system. To find the parameters involved in the 6 state Markovian model for I Ca,L , a non-linear least-squares method  was used to obtain the solution of the system of non-linear ordinary differential equations. Specifically, we have employed an algorithm given by Lau . The numerical integration scheme used to solve the full set of forty two 1st-order differential equations describing the dynamic model was the Merson-modified Runge-Kutta 4th-order method [20, 21] with a conservative fixed time step, chosen small enough to allow the local truncation error to be of fourth order. The explicit finite difference scheme was used to numerically solve the Laplacian equations of Ca2+ diffusion in the cleft space. Detailed numerical methods are similar to those presented by Smith et al. . The results were visualized using Matlab by Mathworks and Origin by Microcal Software.
Our objective was to develop a model of the rat ventricular cell which could be used to explain Ca2+ signaling at the nanoscale level of the dyad and integrate the contributions of many dyads to produce a Ca2+ transient and continuous Ca2+ balance at the whole-cell level. Therefore, we start with a broad discussion of the elements of the DCU and its Ca2+ supply (the jSR), and continue with a progressively more detailed description of the whole cell model. It is important to note that all Ca2+ concentrations discussed in the model pertain to unbound Ca2+ unless specified.
We assume that a continuous membrane barrier exists between the cytoplasm and the external bathing medium (Figure 1A; Figure 2), which consists of two components: a surface sarcolemma (SL) (M FreeSL in Figure 3) free of any sub-membrane contact with the junctional sarcoplasmic reticulum (jSR) and the remnant membrane (M JunctionalSL in Figure 3) that does make contact with the jSR via a dyadic space (nanodomain) (Figure 1A). These membrane components have the same basic plasma membrane, but differ in content with regard to total membrane surface area, type and distribution of transmembrane ion channels, ATPase pumps and exchangers, as well as their functional coupling with a dyadic space. Ultrastructural information from several cardiac preparations including the rat ventricular cell has been compiled by Bers , which can be used to estimate the percentage of the cell membrane in contact with a dyadic space, for either the free surface plasmalemma or for the transverse tubule (TT) which brings the extracellular medium to the plasma membrane of the dyadic coupling unit. Thus, the bounding membrane is divided into two lumped parts (free and coupled) based on the existence of sub-membrane coupling to a dyadic space (Table 1). A portion of membrane could be part of a transverse tubular membrane, but if there is no dyadic coupling involved, that membrane would be classified as belonging to the free surface plasmalemma. Another portion of membrane might be part of the bounding outer surface of the cylindrical cell and yet have submembrane coupling to a dyadic space. In this case, it would be classified as belonging to the coupled category. Table 2 gives values for volumes of the fluid compartments (shown in Figure 3) assumed for the rat ventricular cell, which are largely based on measured data from rat ventricular myocytes .
Channel and Exchanger Distribution
Recent research has also shown that besides L-type Ca2+ channels, Na+/Ca+ exchanger activity is also found predominantly in the T-tubules of rat ventricular myocytes . Our model configuration reflects this finding in that the tubular fraction of I Ca,L ,I NaCa and I NaCs channels facing a unitary dyadic space are denoted as i Ca,L,TT (source of trigger Ca2+ into a unitary dyad), i NaCa,TT and i NaCs,TT respectively (Figure 2). The free sarcolemmal component of these same channels are denoted as I Ca,L,SL , I NaCa,SL and I NaCs,SL respectively. I Ca,L,TT is the total current entering through L-type Ca2+ channels via all the dyadic units (N dyad ). We define the total L-type current I Ca,L as the combination of I Ca,L,TT and I Ca,L,SL (i.e., I Ca,L = I Ca,L,TT + I Ca,L,SL ). I Ca,L in our model is mostly (90%) from the L-type Ca2+ current in the T-tubules, since Kawai et al. , found L-type current to be highly concentrated (9-fold) in the T-tubules (I Ca,L,TT ) vs. the cell surface sarcolemma (I Ca,L,SL ) of rat ventricular myocytes. We described the I Ca,L channel using a 6-state Markovian model as shown in Figure 4A. The distribution of I NaCa and I NaCs correspond to that of I Ca,L channel in order to ensure Ca2+ and Na+ ion balance.
Since our study is focused on voltage clamp testing of Ca2+ transients in rat ventricular cells, we assume that the majority of Na+ and K+ channels are blocked by either the holding potential used (-40 mV) or appropriate blocking agents. Thus, these channels are not modeled, and we assume that only the dihy-dropyridine (DHP) -sensitive Ca2+ channels, the electrogenic pumps, Na+/Ca+ exchangers and Na+/Cs+ pumps expressed in the free and/or coupled SL membranes contribute to the voltage clamp response. Table 3 provides values for various parameters used to model the ion transport across the sarcolemmal membrane.
The SR Fluid Compartment
The SR is an intracellular organelle that consists of two lumped fluid compartments (the jSR and LSR) that communicate (Figure 1A; Figure 3). Like the sarcolemma, the bounding membranes of the jSR and LSR are differentiated regarding their ionic current content and degree of coupling with the sarcolemma. With regard to ionic currents, the LSR membrane has a thapsigargin-sensitive SERCA pump for pumping Ca2+ into the LSR lumen against a concentration gradient. In contrast, the jSR membrane contains an outwardly directed ryanodine (Ry)-sensitive channel for Ca2+ release from the jSR to the dyadic space. The jSR fluid compartment contains the Ca2+ binding protein calsequestin as well as the proteins triadin and junctin, which interact with the ryanodine receptor (RyR) and calsequestrin. This co-located configuration of the RyR receptor, along with the proteins calsequestrin, triadin and junctin which exist on the luminal side of the jSR membrane, constitutes a jSR Ca2+ release regulating mechanism called the luminal sensor (Figure 4B; Figure 5). The protein-protein interaction between them plays an important role in regulating the open-state of the RyR Ca2+ release channel . A six-state Markovian scheme (Appendix A1, Equations 87-92) is used to describe the dynamics of this interaction and it is called the SR luminal Ca2+ sensor. Figure 4B shows a functional diagram of the luminal sensor and its output state is shown connected to the four-state RyR model. Specifically, the sensor adjusts Ca2+ dependent rate functions within the ryanodine receptor model, which affects the open probability P o of the SR Ca2+ release channel.
With regard to coupling, the DHP and Ry-sensitive Ca2+ channels are assumed to be located on opposite sides of the small dyadic fluid space (nanodomain) as shown in Figure 6, and coupled functionally by a CICR mechanism. The dyadic space is assumed to be in fluid communication with the cell cytoplasm via a restricted diffusion region. In contrast, the LSR is not functionally coupled to the sarcolemma, but rather is in contact with the cytoplasm via the SERCA pump (as shown in Figure 3). Table 4 provides values for parameters used to model the intracellular ion transport.
The Dyadic Coupling Unit (DCU)
In describing this functional unit it is necessary to provide progressively more detailed descriptions of the component elements of the individual dyad, particularly the geometrically opposed DHP and Ry-sensitive Ca2+ channels, as well as the geometry and buffering properties of the small dyadic space. As will be shown, we have described the dynamics of both types of Ca2+ channels using Markovian state models (Figure 4) which include features such as Ca2+ mediated channel inactivation, a graded CICR process with a "calcium gain" of approximately 6-7, and two-dimensional Ca2+ diffusion within the dyadic space. Crank  discusses diffusion problems in a two-phase heterogeneous medium and shows that diffusion through a system of barriers (RyR feet structures in the dyadic cleft space) can be approximated by diffusion in the same region without barriers but with a reduced effective diffusion coefficient. We hence take this approach in modeling the Ca2+ diffusion by solving the 2-D Laplacian equation Appendix A3 (Equations 147-150) in the DCU without explicitly accounting for local potential fields. The DHP-sensitive i Ca,L,TT channel brings in trigger Ca2+ (0.1 pA which is of the same order as measured by Wang et. al. ) causing a sparklet (a local increase in Ca2+ concentration at the mouth of the channel). This trigger Ca2+ causes a release from a cluster of opposing RyR channels, causing a spark. This combined release from a cluster of RyR channels causing a spark is represented as the release from a unitary RyR channel (i RyR ) in our model (shown in Figure 6). The characteristics of elemental Ca2+ release from a unitary RyR channel in our model agrees with data in terms of amplitude which is of the order of 3 pA (reported by Cheng et al.  and Blatter et al. ) and duration (full duration at half maximum (FDHM)) which is of the order of 50 ms (reported by Zima et al. ).
The single DCU in our model represents the lumped activity of a large number of individual dyads (e.g. 10,000), and it is charged with the task of forming the cytosolic Ca2+ transient (hence mechanical contraction) each beat of the cardiac muscle cell. In response to tonic application of voltage clamp pulses, the DCU strongly depends on an adequate supply of Ca2+ from the SR. The measurements of Diaz et al.  show that, although trigger current may be supplied regularly by tonic voltage clamp pulses, there is an inherent steady-state dependence of the magnitude of Ca2+ release on the particular value of SR Ca2+ content (i.e., there is a relationship between SR Ca2+ content and peak [Ca2+] myo ; Figure 4 Diaz et al. ). This plot gives us a glimpse of the input-output relationship of the dyadic coupling unit and indicates that SR Ca2+ content is an important controlling variable for the CICR process implemented by the DCU model.
L-Type Ca2+ Current
A multiple state characterization of I Ca,L in rat ventricular myocytes has been reported previously by our group . The gating scheme used in this I Ca,L model has an additional high voltage state C 6 dhpr as shown in Figure 4A which is introduced to reproduce I Ca,L tail currents. Upon voltage-dependent activation, the channel achieves the primary open state O 2 dhpr . The degree of opening is enhanced in the presence of activated calmodulin-dependent kinase (CaMKII act ) [31–33] which is known as Ca2+-dependent facilitation (CDF). CaMKII has been shown to tether to the I Ca,L channel  functioning as a Ca2+ signalling sensor for facilitation. The open probability of the I Ca,L channel is also increased in the presence of activated calcineurin (CaN act ) . These two effects are modeled as shown in Figure 4A, where is a function of both CaMKII act and CaN act besides voltage. The gating scheme features two different pathways for inactivation of the open state. The pathway O 2 dhpr ↔ C 5 dhpr accounts for voltage-dependent inactivation, whereas the pathway O 2 dhpr ↔ C 4 dhpr ↔ C 5 dhpr with Ca2+-calmodulin (Ca 4 CaM) dependent rate constant accounts for the fast and slow phases of Ca2+-dependent inactivation (CDI) . Transitions to the closed state via both the inactivation pathways are suppressed in the presence of CaMKII act and CaN act to allow CDF. A 2-state Markovian model allows Ca2+ mediated interaction between calmodulin and the I Ca,L channel. As shown in Figure 4A, state S 2 dhpr , which denotes Ca 4 CaM bound to the IQ-motif of the I Ca,L channel, modulates calmodulin dependent Ca2+-induced inactivation. This is in agreement with the findings of Nikolai M. Soldatov . Our approach to modeling the effects of the proteins CaM, CaMKII, and CaN on the I Ca,L channel was based on the premise that they are co-localized with the channel itself . Most of the beat-to-beat modulation is produced by CaM and CaMKII, whereas CaN is constitutively active in the dyad . Tables 5, 6 and 7 provide values for the rate constants used to model Ca2+/CaM binding, CaM buffering, Ca2+/CaM/CaMKII as well as Ca2+/CaM/CaN interactions.
Our previous study  was focused on the characterization of the I Ca,L channel under conditions of low [Ca2+] in the dyadic space and myoplasm. In fact, Ca2+- release from the jSR was blocked by administering a relatively high dose of ryanodine (20 μ M) in all experiments. Therefore, to study the additional influence of SR Ca2+ release on I Ca,L , we modified the original I Ca,L model to better characterize the process of Ca2+-dependent inactivation. In the modified I Ca,L model, majority of the structure for the 6-state dynamic scheme remains the same, however changes have been made to the voltage-dependent inactivation rate function (), and the Ca4CaM dependent inactivation rate function (). The specific formulas for the Markovian state equations and the modified and functions used in this study are given in Appendix A1 (Equations 4-50). Our adjustments consist of reducing the contribution from voltage-dependent inactivation process and strengthening the Ca4CaM-dependent inactivation process. The rate constants and are constructed in order to provide a re-excitation window during the return of the clamp voltage to the resting potential. With these adjustments, the Markovian state description for I Ca,L can provide good fits to measured I Ca,L data under both test conditions (presence and absence of ryanodine -sensitive Ca2+ release) as well as produce tail currents during repolarization from large clamp voltages (≥ 40 mv).
Calcium Buffering in the Dyadic Space
Previous modeling work by Post et al. , and Post and Langer  considered the effect of Ca2+ binding sites on the inner sarcolemmal leaflet. Following these authors, we included low-affinity (k d = 1.1 mM) and high-affinity (k d = 1.3 μ M) Ca2+ binding sites on SL wall boundary of cylindrical dyadic space. The presence of these membrane bound sarcolemmal Ca2+ binding sites in our model has significant physiological implications, in that it prevents local dyadic Ca2+ concentration near the "mouth" of DHP-sensitive I Ca,L channel from becoming excessively high. Allowing such a condition to occur can cause a reversal of the I Ca,L current, which does not occur physically during normal jSR release. Addition of these SL Ca2+ binding sites does not significantly slow the build-up of Ca2+ within the dyadic space, nor is the Ca2+-induced Ca2+ release mechanism affected .
Ca2+ Release Channel
The gating characteristics of the Ry-sensitive release channel are not only modulated by the dyadic Ca2+ concentration at its mouth but also the jSR Ca2+ concentration via the luminal sensor. Several RyR gating schemes have been deduced from isolated RyR currents measured in lipid bilayers, including: (1) the 4-state scheme developed by Keizer and Levine  with 4 Ca2+ binding for activation and 3 Ca2+ binding for inactivation to explain the "adaptation" of the RyR observed by Gyorke and Fill ; (2) the 6-state scheme suggested by Zahradnikova and Zahradnik , which allowed opening the RyR channel upon binding of a single calcium ion; and (3) the Markov model proposed by Keizer and Smith , which can be dynamically switched among the six, five and four-state representations during the simulation as Ca2+ levels vary. Stern et al.  demonstrated that none of these schemes yielded stable local control of SR release, even with extensive adjustment of parameters. Stern et al.  further reported that all schemes resulted in models that manifested local instability, as indicated by failure of release to terminate after activation, or global instability caused by spontaneous activation by resting [Ca2+] myo . Since many of the kinetic gating schemes derived from lipid bi-layer data fail to support stable E-C coupling in simulations, he concluded that the RyR gating process in situ may differ considerably from that in bi-layers.
Our gating scheme is patterned after the release channel used in the model of Stern, Pizarro and Rios  (Figure 4B), where the channel is assumed to have four states: rest (closed), activated (open), inactivated (closed) and refractory (closed). The activation gate is opened by the simultaneous, cooperative binding of two Ca2+ ions, whereas inactivation depends on the binding of a single Ca2+ ion. Four additional features have been built into our model: (a) the crucial role of the luminal SR Ca2+ sensor (Figure 5) in assisting the inactivation of the RyR channel is modeled via the dependence of all the rate constants in the 4-state RyR model, on the degree of interaction between the RyR and the proteins triadin and junctin (Figure 4B); (b) CaMKII act dependent enhancement in RyR release [47–50] is modeled via the rate functions , , and ; (c) a stronger Ca2+ dependent inactivation of the RyR channel (at a fixed depolarization level), is adopted to reflect recent observations that the inactivation of the RyR depends on the high local Ca2+ concentration consequential to their own Ca2+ release ; and (d) inactivation of the RyR at different depolarization levels is made dependent on local Ca2+ concentration (i.e., [Ca2+] ryr at the "mouth" of the RyR channel on the dyadic side) as per Wier et al.  and Zucchi et al. . Our initial studies indicated that the repriming rate () in the RyR gating scheme of Stern et al.  was quite large, which can lead to a saturated open probability of the Ca2+ release channel. This occurs during the later phase of channel inactivation, where the large value of tends to reactivate the channel, resulting in saturated Ca2+ release. Therefore, we utilized a value of that is 10% of that used by Stern et al. , and further assumed that the unitary permeation flux of the jSR release channel is proportional to jSR luminal Ca2+ concentration ([Ca2+] jSR ). The specific equations for the Ca2+ release model are given in Appendix A1 (Equations 73-92).
Luminal RyR sensor
The RyR Ca2+ release is modulated by a multi-molecular Ca2+ signalling complex which is localized to the junctional SR [54–56]. This complex consists of the ryanodine receptor (RyR) which functions as a Ca2+ conducting pore [57, 58], calsequestrin (CS) which acts as the Ca2+ binding protein [59, 60], and the junctional SR transmembrane proteins triadin  and junctin . It was known previously that the proteins triadin and junctin anchor calsequestrin to the ryanodine sensitive receptor . More recently, it has been observed that the protein-protein interactions between triadin, calsequestrin and RyR modulate sarcoplasmic reticulum calcium-release in cardiac myocytes . Triadin and junctin are structurally homologous proteins , but the functional differences in their roles are unclear at present. Therefore we have refrained from modeling these proteins separately with unique roles, and will hereafter only mention triadin. The RyR model proposed by Shannon et al.  takes a heuristic approach towards free SR Ca2+ concentration dependent luminal control of the RyR channel. Our detailed biophysical model of the luminal sensor is strongly based on the recent findings of Terentyev et al.  which uncovers complex [Ca2+] jSR ) dependent, CS mediated mechanistic interaction of the protein triadin with the RyR channel.
The triadin protein facilitates SR Ca2+ release by sensitizing the RyR to activation by the trigger current I Ca,L,TT . This is incorporated in our model by allowing the rate constants in the 4-state Markovian model for the RyR (Figure 4B-ii) channel to be functions of the activated state A 1 ls , which represents the degree of binding between triadin and RyR, in the 6-state model for the luminal sensor (Figure 4B-i).The degree of triadin assisted anchoring of CS to the RyR channel [61, 62] is denoted by the state I 2 ls . Triadin is also known to exist in a form bound to CS alone  denoted by the state I 3 ls and in its free unbound form represented as I 4 ls in the model. CS which is known to exist in a Ca2+ bound form modeled by B 6 ls also modulates SR Ca2+ release by influencing the open probability of the RyR channel [65–67], via interactions with Triadin [54, 66, 25]. The degree of unbound CS is denoted by B 5 ls in the model.
During the diastolic period the relatively large concentration of available free Ca2+ in the SR results in most of the CS being bound to Ca2+, decreasing the degree of interaction between CS and triadin. This enables a strong interaction between the available unbound triadin and the RyR channel increasing its propensity for trigger Ca2+ induced SR release. Our model incorporates this property by facilitating movement of states towards A 1 ls and B 6 ls in the presence of large SR Ca2+ concentration. Following SR Ca2+ release, a reduced luminal Ca2+ concentration in the jSR causes an increase in the amount of free calsequestrin (B 5 ls in the model) available to bind with triadin. This results in a decrease in the extent of interaction between triadin and RyR (A 1 ls in the model), thus inhibiting the RyR channel and leading to robust termination of SR Ca2+ release in cardiac myocytes . This release termination mechanism is incorporated in our combined RyR-Luminal sensor model (Figure 4B and Appendix A1, Equations 73-92).
Due to the lack of in-vivo measurements on individual state transitions, the rate constants in our novel luminal sensor model (Table 8) are chosen such that (i) the triadin mediated sensitization of the RyR channel (via A 1 ls ) provides adequate peak RyR release which translates into the upstroke velocity of the cytosolic Ca2+ transient; (ii) the luminal sensor mediated RyR channel inactivation (via A 1 ls ) causes timely release termination resulting in a cytosolic Ca2+ transient duration, physiological for a rat ventricular myocyte; (iii) the rate of post-release RyR recovery results in appropriate channel refractory characteristics. In the case of channels where there is an interaction between ion flow (which is not at equilibrium) and its gating mechanism the microscopic reversibility criteria does not hold true [68, 69]. The RyR channel which is solely modulated by [Ca2+] (both on the dyadic side ([Ca2+] ryr ) as well as the luminal side ([Ca2+] jSR )) experiences a strong interaction of Ca2+ flow through itself and its gating mechanism described as [Ca2+] ryr induced self-inhibition and the luminal sensor induced inactivation. Hence, the rate constants of the luminal sensor model (Figure 4B-i) are not constrained by microscopic reversibility criteria. However, a stability constraint in the form of, the sum of probabilities of all possible states corresponding to triadin (A 1 ls , I 2 ls , I 3 ls , I 4 ls ) and CS (I 2 ls , I 3 ls , B 5 ls , B 6 ls ) being equal to one is explicitly imposed Appendix A1, Equations 87-88) on the model of the luminal sensor in order to avoid run-off.
Ca2+ Buffering in Myoplasm and SR
Ca2+ buffers play an important role in sequestering a fraction of the total Ca2+ released during E-C coupling and contraction. These buffers include: (a) calmodulin (CaM), which is assumed to be uniformly distributed in the myoplasm and dyadic space; (b) troponin in the bulk myoplasm; and (c) calsequestrin (CS) in the jSR. The dyadic space has been shown to be accessible to calmodulin (CaM), but mostly inaccessible to fluorescent dyes [70, 71]. Therefore, we consider the dyadic space filled with calmodulin, but not fluo-3. This provides a more direct calcium communication pathway between the DHP and Ry receptors. The rate constants for Ca2+ binding to calmodulin were based on a model from , whereas those for Ca2+ binding to troponin were taken from Potter and Zott . Rate constants used to describe Ca2+ binding to calsequestrin were based on the study of Cannell and Allen , whereas those for Ca2+ binding to troponin-Mg complex were adopted from Lindblad et al. .
It is well recognized that fluorescent indicator dyes introduced into the cytosol also act as Ca2+ buffers , even at submillimolar concentrations . In our simulations, we have used 100μ M fluo-3. We assume that the quantity observed experimentally as a "calcium concentration" signal is actually the calcium complexed with fluo-3, or [CaF3]. The differential equation describing the change in [CaF3] with time is given in Appendix A2 (Eq. 138) which follows from Shannon et al. .
Cytosolic Ca2+ is pumped into the LSR (Figure 1A) by a Ca2+-ATPase, which lowers [Ca2+] myo and helps to induce relaxation in cardiac muscle. The transport reaction involves two Ca2+ ions and one ATP molecule , and it is represented by the description given for I cyt,serca and I serca,sr in Appendix A1 (Equations 51-63). The model used for the uptake pump is adopted from Koivumaki et al.  and takes into account both the forward flux of Ca2+ from the cytosol to the LSR lumen and the backward flux from the LSR to the cytosol along with the Ca2+ buffering action of the SERCA protein. The phospholamban (PLB) to SERCA ratio has been fixed to 1.0, assuming almost equal availability of both the proteins. CaMKII act affects the SERCA pump via direct phosphorylation assisting in enhancement of SR Ca2+ transport by increasing the pumping rate  and indirectly via phosphorylation of PLB  relieving the inhibition caused by PLB on the SERCA pump in turn increasing the sensitivity of the pump for Ca2+ uptake. These two effects are modeled by allowing the rate constants for Ca2++ binding to/release from the SERCA pump as well as the rate constant for phosphorylation of PLB to be a function of CaMKII act in the cytosol. The activating role of CaN in modulating the SERCA pump [80, 81] is accounted for in our model by allowing the rate constant for phosphorylation of PLB to be dependent on available CaN act in the cytosol. The relative regulatory role of the enzyme protein kinase A(PKA) is fixed at 0.1, as its modulatory effect is beyond the scope of this study. The current I cyt,serca dictates the transport of Ca2+ between the cytosol and the SERCA protein. Similarly, the current I serca,sr dictates the transport of Ca2+ between the SERCA protein and the LSR. The difference in these Ca2+ currents accounts for the Ca2+ buffered by the SERCA protein. The jSR is subsequently refilled by Ca2+ diffusion from the LSR. The differential equations describing the Ca2+ balance and particularly the transfer between two separate SR compartments (jSR and LSR) are provided in Appendix A1 (Eq. 72).
Ca2+-Extrusion via Sarcolemmal Ca2+ Pump
Although the sarcolemmal Ca2+-pump has a high affinity for [Ca2+] myo , its transport rate is far too slow for it to be an important factor in Ca2+ fluxes during the cardiac cycle. It might, however be more important in long-term extrusion of Ca2+ by the cell. Our model of the plasma membrane Ca2+ ATPase pump current is adopted from Sun et al. . We have used a constant value of half activation constant k mpca = 0.5μ M in our model on the basis of measurements by Caroni et al .
Ca2+-Extrusion via Na+/Ca2+ Exchanger
In mammalian cardiac cells, it is generally accepted that the Na+/Ca2+ exchanger has a stoichiometry of 3Na+:1Ca2+ . I NaCa (again, the combination of I NaCa,TT and I NaCa,SL ) is important in removing Ca2+ during twitch relaxation, in competition with I up . A simple thermodynamic Na+/Ca2+ exchanger current model  may be sufficient to predict the direction of Ca2+ transport by Na+/Ca2+ exchange and the driving force, however the amplitude is subject to kinetic limitations (depending on substrate concentrations). A more comprehensive Na+/Ca2+ exchanger current equation  is adopted in our model, which includes factors for allosteric Ca2+ activation and the transport for Na+ and Ca2+ inside and out. The maximal flux through the exchanger V max is estimated to ensure that the Ca2+ ion transport (which is voltage dependent) via the Na+/Ca2+ exchanger matches the influx of Ca2+ via the I Ca,L channel [86, 87], maintaining whole cell Ca2+ homeostasis.
The Na+/K+ pump, helps in maintaining homeostasis of the intracellular Na+ ion concentration. ATPase activity powers the pump, as it generates an outward Na+ flux and an inward K+ flux with a Na+ to K+ stoichiometry of 3 to 2 . However, in our experimental protocol, external solution in the bath was normal Tyrode (1 mM Ca2+) with Cs+ substituted for K+ in order to block the inward rectifier K+ current. The internal solution in the pipette contained Cesium aspartate supplemented with 20 mM CsCl, 3 mM Na2ATP , 3.5 mM MgCl2 and 5 mM HEPES.
Activation of the electrogenic sodium pump in mammalian non-myelinated nerve fibres , skeletal muscle  and rat brain cells  by Cs+ is reported in the literature. It is observed that the late effects of reducing extracellular K+ concentration ([K+] o ) to 0 mM in mammalian cardiac muscle can be prevented by including appropriate concentrations of other activator cations of the Na+/K+ pump such as Cs+ in the 0 mM [K+] o bathing solution . Monovalent cations (including Cs+) were also added to K+ free bathing solution to reactivate the sodium pump in guinea-pig ventricular myocytes . The effectiveness of Cs+ as an external cation in activating the electrogenic sodium pump is known to be lesser than potassium .
We have represented I NaCs in our model by the expression for Na+/K+ pump formulated by Linblad et al.  replacing K+ ion concentrations with the Cs+ ion concentration. While ensuring whole cell Na+ ion balance, the peak Na+/Cs+ pump current is modified to be one-sixth to account for the decreased potency of the cation Cs+ in activating the pump. The voltage-dependence of I NaCs is adopted from the data on Na+/K+ pump from Hansen et al. .
The DCU is a fundamental element in the mechanism of CICR. The sequence of events resulting in CICR is triggered by the Ca2+ entering through the I Ca,L channel. The characteristics of this channel are hence examined in detail to understand its voltage and Ca2+ dependent behavior. The ability of the I Ca,L channel to facilitate graded RyR release is noted. The high gain associated with RyR release is quantified. This is followed by a study of the properties of the RyR channel with a particular emphasis on its inactivation mechanism. The Ca2+ mediated interaction between these channels is investigated in detail. This is followed by an effort to understand the role of [Ca2+] jSR in CICR. In particular, the relationship of the peak cytosolic calcium transient and the SR Ca2+ content is examined. The effect of modulatory agents like caffeine and thapsigargin on SR release/uptake is studied to understand the significance of these mechanisms in facilitating a normal CICR.
L-type Ca2+ current (ICa,L)
The L-type DHP-sensitive Ca2+ channel has a key role in initiating CICR. Hence a thorough analysis of its activation and inactivation mechanisms is considered here. Though the activation of the I Ca,L channel is solely voltage dependent once activated, the inactivation of the channel is influenced not only by the trans-membrane voltage but also the Ca2+ concentration in the vicinity of the channel . The relative con-tribution of the voltage-dependent and Ca2+-dependent inactivation pathways in the I Ca,L channel model can be studied by selectively blocking each pathway in the model. Figure 7 shows the model-generated waveform for I Ca,L under three conditions; (i) control case with both voltage and calcium dependent in-activation pathways intact; (ii) ryanodine applied to allow only the voltage-dependent inactivation and the Ca2+-dependent inactivation via I Ca,L self-inhibition in the absence of RyR release; and (iii) Ca2+ substituted with Ba2+ to facilitate only the voltage-dependent inactivation pathway and block all Ca2+ dependent inactivation pathways (k 24 dhpr set to zero). The protocol used in case (i) was a 50 ms pulse with clamp voltages ranging from 10 mv to 40 mv in steps of 10 mv from a holding potential of -40 mV. In case (ii) and (iii) the pulse duration was increased to 200 ms to replicate the data.
Measured data obtained using the same voltage clamp pulse under all three conditions (Control, Ryan-odine application and Barium substitution) is shown for comparison (data obtained from Dr. Palade's lab was pre-processed to eliminate transients produced by changing clamp voltage in order to obtain model fits using a non-linear least-squares method ). The degree of I Ca,L channel opening is known to be enhanced in the presence of CaMKII act  which is known as Ca2+-dependent facilitation (CDF). According to Bers et al.  the positive regulation of the channel by CaMKII act requires Ca2+ influx (it is not seen when Ba2+ is the charge carrier as in Figure 7D and is more strongly apparent when local Ca2+ influx is amplified by SR Ca2+ release as in Figure 7A). Increase in CaMKII act seems to play an important function in enhancing I Ca,L peak amplitude , suggesting a critical role for Ca2+-dependent facilitation. Our model generated results also indicate that a stronger CDF (in Figure 7A) in the presence of RyR Ca2+ release causes an enhancement in I Ca,L peak amplitude (compare the model fits to data in Figure 7A (presence of RyR release) with Figure 7B (no RyR release) and Figure 7C (Ca2+ substituted by Ba2+)). It is important to note that the peak of the I Ca,L current is reached after the RyR open probability reaches its maximum value owing to the faster dynamics of the RyR channel. Our simulations indicate that following the inward current peak, Ca2+-dependent inactivation (CDI) is much faster and dominates the response for the duration of voltage clamp (30 ms <t <80 ms). In addition, a comparison of Figs. 7A and 7B shows that, Ca2+-dependent inactivation (CDI) caused by Ca2+ release from Ry-sensitive channels is much more significant than the self-inhibition produced by Ca2+ influx via I Ca,L itself. The graded behavior of the cytosolic Ca2+ transient is shown in Figure 7C. Voltage dependent inactivation (VDI) is relatively slow compared with CDI and is best seen in Figure 7D where all Ca2+ inactivation effects have been blocked by Ba2+ substitution for Ca2+. Under RyR blockade (Figure 7B), the relatively slow VDI has its major effect during the late phase of the long voltage clamp pulse (e.g., beyond 100 ms) and in a time range where CDI is relatively constant. Quantitative analysis of movement of states (during the depolarizing pulse duration of 50ms) via different pathways in the six state Markovian model shows that 68.24% of the total inactivation of I Ca,L,TT is via the calcium dependent O 2 dhpr -C 4 dhpr pathway, owing to the large Ca2+ concentration which the channel is exposed to in the dyad. In contrast, only 19.21% of the total inactivation of I Ca,L,SL is via the calcium dependent O 2 dhpr -C 4 dhpr pathway, because of the low cytoplasmic Ca2+ concentration in the vicinity of the channel. Considering the total I Ca,L current (I Ca,L,TT + I Ca,L,SL ), around 63.34% of the I Ca,L channel inactivation is via the Ca2+ dependent pathway.
Figure 8A shows a comparison of model generated and experimentally obtained I Ca,L data (obtained from Dr. Palade's lab and pre-processed to eliminate transients produced by changing clamp voltage in order to obtain model fits using a non-linear least-squares method ) from a rat ventricular myocyte (different from the cell used to obtain data in Figure 7) at negative (-30 mv ≤ v ≤ 0 mv) clamp voltages. The corresponding graded behavior of the cytosolic Ca2+ transient is shown in Figure 8C. Figure 8B shows the well known  bell-shaped dependence of the peak [Ca2+] myo on the clamp voltage.
Figure 9 shows model generated normalized I Ca,L obtained by a voltage clamp to 10 mv from a holding potential of -40 mv using the following protocols. Case 1: Normal release is allowed where, following CICR, a strong Ca2+ dependent inactivation inhibits the I Ca,L current as seen in the plot numbered 1 (data obtained from Dr. Palade's lab). Case 2: Ryanodine applied to block RyR release; Ca2+ entering via the I Ca,L channel causes Ca2+ induced inactivation, although the magnitude of inhibition is far less than in the control case. Case 3: Ba2+ substitution for Ca2+ to completely suppress the Ca2+ inactivation mechanism; the only inactivation pathway present is the slow voltage dependent pathway , which causes substantially reduced recovery compared to cases 1 and 2. Besides the inactivation, it can be seen from Figure 9 that the rate of channel activation, which is evident in the slope of the individual plots during the initial activation phase, is faster in the presence of RyR Ca2+ release (a result of enhanced CDF in the presence of elevated Ca2+ levels). This highlights the important role of Ca2+ in regulating I Ca,L channel opening and thus controlling the amount of extracellular Ca2+ entering the cell.
ICa,L,TT - dependent Graded SR-release
One of the most important characteristics of the CICR mechanism is the modulation of the RyR Ca2+ release based on the amount of trigger calcium delivered via the DHPR channel. The voltage-controlled I Ca,L channel is the link between the extracellular excitation and the intracellular Ca2+ release. As shown in Figure 10A, our model reproduces graded release. It is important to note that the onset of the Ca2+ transient is also modulated based on the magnitude of the trigger calcium available to initiate release, as shown in Figure 10A-iii. This is due to the fact that the rate of increase in the open probability of the Ry-sensitive Ca2+ release channel is controlled by the amount of trigger Ca2+ present at the 'mouth' of this channel, which in-turn is graded by the amount of trigger Ca2+ entering the DCU by means of the DHP-sensitive I Ca,L,TT channel. This mechanism is incorporated in our model by allowing the rate constants in the 4-state Markovian model of the Ry-sensitive Ca2+ release channel (Figure 4) to be [Ca2+] ryr dependent. Small changes in i Ca,L,TT cause modulation of pre-release [Ca2+] ryr , which dictates the propensity of the RyR channel for a Ca2+ release. The pre-release values of the rate constants in the 4-state Markovian model of the RyR channel (nonlinear dependence on [Ca2+] ryr as shown in Appendix A1, Equations 79-84) along with the diastolic [Ca2+] jSR (which is kept constant in Figure 10) set the peak RyR open probability achieved by the RyR channel as well as the peak [Ca2+] myo . The use-dependent adaptation of the RyR channel  is reflected in its non-linear response to trigger Ca2+. As shown in Figures 10B-iii, iv, the delay in the onset of cytosolic Ca2+ transient closely follows the modulation of the onset of RyR release. Though occurring in separate compartments, the peak of cytosolic Ca2+ transient also tracks the maximum value attained by the open probability of the Ry-sensitive Ca2+ release channel. Here, the clamp voltage was held constant at 10 mv to avoid the interference of voltage dependent Ca2+ transport via the Na+/Ca2+ exchanger which is co-located  with the I Ca,L channel in the dyad and the scaling down of I Ca,L,TT corresponds to a fast flicker block of the channels by dihydropyridine.
The slow rising foot that precedes the rapid upstroke (trace at position 0 in Figure 10B-i and B-ii) is the contribution from a single sparklet (Figure 6), which is the result of Ca2+ release from a single i Ca,L,TT channel in our model (0.1 pA), which is of the same order as reported by Wang et. al.  bringing trigger Ca2+ into a unitary dyad. This L-type Ca2+ channel (LCC) triggers release from a cluster of RyR channels causing a Ca2+ spark (Figure 6), which results in the rapid increase in [Ca2+] dyad that follows the foot (shown in Figure 10B-i and B-ii). As shown in Figure 10B, the latency from the onset of the sparklet foot to the triggered Ca2+ spark increases with decrease in the magnitude of i Ca,L,TT . The SR Ca2+ release flux underlying a typical Ca2+ spark corresponds to approximately 2 pA [27, 28, 13]. From RyR single channel conductance measurements in lipid bilayer studies , a Ca2+ spark translates into a release from around 4 RyR channels (also reported by Blatter et. al. ). This single Ca2+ spark corresponds to a unitary RyR release (2 pA) in our model. Our model does not attempt to reproduce the stochastic kinetics of the single channel LCC-RyR coupling; however it mimics accurately the average behavior of this stochastic process which is the net Ca2+ release flux into the cytosol causing the Ca2+ transient. This on/off stochastic nature of the coupling was used earlier to explain the release termination via stochastic attrition. However, it is now known  that the luminal sensor plays a fundamental role in an active extinguishing mechanism  that effects a robust [Ca2+] jSR - dependent closure of the RyR channel. This mechanism for inactivation of the Ry-sensitive Ca2+ release channel is accounted for in our model.
Cytosolic Ca2+ transient is also graded by the duration of the I Ca,L,TT trigger current controlled by the voltage clamp pulse duration. Our model also shows that triggered release can be prematurely stopped by rapid repolarization [100–102]. The effect of depolarization duration on the time course of the cytosolic Ca2+ transient is indicated in Figure 11 where the duration of the voltage clamp pulse is decreased from 80 ms to 5 ms (while keeping the clamp voltage constant at 10 mv) resulting in premature stoppage of release and thus a decrease in peak cytosolic Ca2+ transient. This effect is a combined result of (i) the modulation of release due to pulse duration dependent change in the amount of trigger Ca2+ entering the dyadic coupling unit (DCU) and (ii) the pulse duration dependent change in the relative role of the Na+/Ca2+ exchanger co-located in the DCU.
High Gain of Ca2+ Release
Besides the graded release, an extremely valuable characteristic of the CICR process is the high gain associ-ated with it. A small amount of Ca2+ entering the DCU via the I Ca,L,TT channel causes a large release of Ca2+ from the sarcoplasmic stores via the Ry-sensitive Ca2+ release channel on the jSR lumen that interacts with the DCU. This high gain is essential in producing physiological cytosolic Ca2+ transients when 10000 of these DCU's operate in tandem. Two different definitions of the gain or amplification factor due to CICR have been adopted in our model simulations, namely the ratio of:
average integrated RyR flux to average integrated DHPR flux ; and
the peak Ca2+ transient in the presence of CICR to the peak calcium transient in its absence, con-tributed by the trigger calcium alone .
Thus, criterion (1) measures calcium gain observed in the dyadic space, whereas criterion (2) measures gain observed in the cytosolic space.
With regard to the dyadic space, Figure 12A shows the Ca2+ flux through DHPR and RyRs, respectively. Figure 12B shows cytosolic calcium transients under conditions of ryanodine blockade and with RyR Ca2+ release. Based on our measurements, there is approximately a 2 ms delay from the onset of the DHPR influx to the initiation of the RyR release flux.
At a clamp voltage of 10 mv, criterion (1) applied to our simulation yields an integrated flux ratio of 6.39, whereas a gain of 5.75 is calculated using criterion (2). By either method, a CICR amplification factor of approximately 6 is calculated, which is similar to that reported by Stern . This result is also consistent with gain calculations from the measured data of Fan and Palade  on rat ventricular cells. They estimated a gain of approximately 7 using comparisons of the rates of rise of Ca2+ transients in the presence and absence of ryanodine. The model generated results are obtained by using a voltage clamp protocol of a 50 ms step pulse to 10 mv from a holding potential of -40 mv.
The inset in Figure 12A shows the voltage dependence of CICR gain formulated using criterion (1). Our model shows a decline in gain as the clamp voltage is increased from -30 mv to 20 mv [52, 103]. However, any further increase in clamp voltage results in a small increase in gain (Figure 1C, Altamirano et al. ). It is important to note that, with increasing clamp voltage, the decreased ability of the Na+/Ca2+ exchanger (which is co-located  in the dyad) to extrude Ca2+ partially compensates for the declining trigger current, in facilitating SR release and hence assists in increasing CICR gain.
RyR Refractory Characteristics
Given the fact that the CICR mechanism has a high gain associated with it, it could be prone to unstable behavior if not for the refractory characteristics of the Ry-sensitive release channel . It is now known  that the luminal sensor plays a fundamental role in an active extinguishing mechanism  effecting a [Ca2+] jSR dependent robust closure of the RyR channel which is accounted for in our model. The RyR release fails to become regenerative due to the role of the luminal sensor which after a release occurs, forces the RyR channel into an absolute refractory state followed by a relative refractory state as shown in Figure 13C. When the RyR channel is in the absolute refractory period, [Ca2+] ryr drops to a level much below what is caused by a sparklet (trigger Ca2+) hence it robustly avoids re-excitation. This extremely critical refractory feature of the channel enabled by the RyR luminal sensor acts as a protective mechanism against premature Ca2+ release in the wake of a secondary excitation that occurs before jSR can be filled back to the control level. The refractory nature of the channel is caused by the [Ca2+] jSR dependent inhibition induced by the protein triadin from the luminal side of the RyR channel. This restraining lock on the RyR channel assists in reloading the jSR.
A dual stimulus protocol (S1-S2) was employed to study RyR refractoriness. The RyR channel was stimulated initially by stimulus trigger current S1, followed after an interval (T2 in Figure 13A) by an identical stimulus current S2. Each of the two stimuli directed into the dyadic coupling unit is the dyadic component (I Ca,L,TT ) of trace 1 (I Ca,L evoked by a 10 mv clamp pulse) in Figure 7A with a peak amplitude of 6.657 pA/pF and a duration (T1 in Figure 13A) of 50 ms. The stimuli S1 & S2 were kept identical to delineate the effects of RyR refractoriness on CICR. The time interval between the two stimuli (T2) was gradually reduced from 250 ms to 50 ms in steps of 25 ms to observe the effects of partial recovery of the RyR channel.
As seen in Figure 13, the decrease in interval between stimuli causes partial recovery of [Ca2+] level in the jSR which manifests in increasing levels of unbound version of the protein calsequestrin, which tends to bind with triadin. As more triadin binds to calsequestrin, there is less left to bind to the luminal side of RyR to activate the channel for CICR. This signaling sequence mediated by the luminal sensor via the interaction of proteins calsequestrin and triadin, allows only a partial SR Ca2+ release based on the degree of SR Ca2+ content recovery. This is an important feature that helps restore the jSR Ca2+ content at the end of every cycle. The RyR channel in our model has an absolute refractory period of 75 ms and a relative refractory period of 250 ms (Figure 13C). These results are consistent with the refractory period measurements by Sobie et al  on rat ventricular myocytes. When the RyR channel is in the absolute refractory period, [Ca2+] RyR drops to a level much below what is caused by a sparklet (trigger Ca2+) and hence robustly avoids re-excitation. However, when the RyR channel is in the relative refractory period, a partial release (Figure 13C) is possible, despite the RyR receptors affinity for release being low.
Figures 14A and 14B show the steady state [Ca2+] myo and [Ca2+] jSR predicted by the model, with and without a functioning luminal sensor, respectively. The stimulation protocol used is a pulse train of amplitude (-40 mv to 10 mv), duration (50 ms) and frequency of 4.0 Hz which is applied for a period of 2.5 sec. The value of the luminal control of the RyR channel in allowing adequate SR filling can be seen in Figure 14 where, in the presence of the luminal sensor the cell maintains a normal cytosolic Ca2+ transient (Figure 14A-i). This is made possible due to sufficient SR filling (Figure 14A-ii) as a result of RyR refractoriness. However, inactivation of the RyR channel is also known to depend on the high local Ca2+ concentration consequential to its own Ca2+ release . The lack of a luminal sensor forces the RyR channel to solely rely on Ca2+ dependent inactivation mechanism. The resulting inadequate RyR inactivation depletes diastolic [Ca2+] jSR level (Figure 14B-ii) causing the cytosolic Ca2+ transient to diminish to new lowered steady state values (shown in Figure 14B-i). The absence of the luminal sensor is modeled by setting the value of the luminal sensing variable ('var' in the model) to a level consistent with it's normal value under 4 Hz, voltage clamp stimulation. This mimics the case where the luminal sensor is insensitive to changes in [Ca2+] jSR and hence non-functional. It is important to note that the decrease (37.07%) in diastolic SR level (as indicated in Figure 14B-ii) is due to the lack of a luminal control on the RyR channel resulting in a reduced rate of RyR recovery which in turn reflects in a compromised SR filling rate. Hence, the presence of a luminal sensor which monitors the SR Ca2+ content, is key to long term Ca2+ stability in the cell.
CICR modulation by the jSR Ca2+ content
From the refractory characteristics of the RyR channel, it is evident that the RyR release is strongly dependent on the jSR Ca2+ content. In fact, the SR Ca2+ release through the RyR channel depends on: (a) the amount of trigger Ca2+ entering by means of the I Ca,L,TT channel (b) the concentration gradient of Ca2+ between the jSR and the mouth of the RyR channel in the dyadic space and (c) the open probability of the RyR channel modulated by the interaction between the luminal sensor and the RyR protein. The SR Ca2+ released also inactivates the trigger current I Ca,L,TT (as shown in Figure 7), hence causing an indirect self-inhibition of RyR release.
Figure 15 shows a phase plot of RyR open probability (O 2 ryr ) versus the Ca2+ concentration at the mouth of the RyR channel in the dyadic space [Ca2+] ryr constructed from model-generated data for a pulse of amplitude -40 mv to 10 mv and a duration of 50 ms. The pre-release diastolic [Ca2+] jSR is 918 μ M. Following excitation by the trigger current during phase A, the RyR channel begins to open, allowing Ca2+ flux into the dyad and elevating [Ca2+] ryr at the mouth of the RyR channel. The rate of opening of the channel is Ca2+ dependent; hence as [Ca2+] ryr increases, RyR open probability increases first slowly and then ever more rapidly. With this onrush of Ca2+, [Ca2+] ryr rapidly equilibrates with the [Ca2+] jSR , reducing and eventually abolishing the concentration gradient across the channel. [Ca2+] ryr soon reaches its maximum value (T1 in Figure 15) and begins to track the decrease in [Ca2+] jSR , there being no further significant concentration gradient existing between [Ca2+] jSR and [Ca2+] ryr at the mouth of the RyR channel. During phase B, the rate of rise in RyR open probability decreases due to decreasing [Ca2+] ryr and the onset of self-inhibition due to the large values of Ca2+ concentration at the mouth forcing the RyR open probability to attain its maximum value (T2 in Figure 15). This is followed by phase C, where the RyR open probability begins to decrease slowly, initiating channel recovery due to decreasing overall [Ca2+] dyad levels as a result of (a) Ca2+ fluxing out of the dyad into the cytosol, (b) lack of a drive from SR release due to drastically diminished Ca2+ gradient at the mouth of the RyR channel and (c) continued Ca2+ induced self-inhibition. Decreasing levels of free Ca2+ in the jSR force a strong RyR inactivation via the luminal sensor (beginning at T4 and continuing throughout phase D; Figure 15), which in turn causes a sharp decline in RyR open probability. The minimum value of [Ca2+] jSR occurs at T5, and the channel is ultimately closed at T6 (Figure 15). The luminal sensor mediated inactivation assists in robust [Ca2+] jSR recovery.
Cytosolic peak [Ca2+] dependence on SR Ca2+ content
Figure 16A shows three characteristic regions in the plot of the peak [Ca2+] myo vs [Ca2+] jSR . In Region I, where the jSR load is small, two things occur: (a) RyR release is reduced, as is the basal intracellular Ca2+ concentration, and (b) Ca2+ dependent inactivation of the I Ca,L,TT channel is reduced due to the lowered dyadic Ca2+ concentration, which allows for greater Ca2+ entry into the dyadic space via the I Ca,L channel. This increase in trigger current results only in a very small increase in Ca2+ release, because of the inhibition of the Ry-sensitive channel caused by the luminal sensor in response to low SR Ca2+ content. Thus, these factors cumulatively result in only a small linear increase in the peak [Ca2+] myo with increasing SR Ca2+ content in Region I. The behavior in Region I corresponds to traces 7-1 in Figure 16B (shows phase plots of O 2 ryr vs [Ca2+] ryr for different steady state [Ca2+] jSR concentrations), where the peak RyR open probability shows very small gradual increase with increasing levels of diastolic pre-release [Ca2+] jSR for a constant trigger Ca2+ input via I Ca,L . The steady state diastolic [Ca2+] jSR levels being low, the maximum RyR open probability values attained are low. Although the peak [Ca2+] ryr does not increase substantially as the diastolic [Ca2+] jSR levels increase from 317μ M to 673μ M, the area enclosed by the [Ca2+] ryr - O 2 ryr loop (which indicates the amount of SR Ca2+ released into the dyad) increases exponentially due to increasing peak RyR open probability combined with the increasing difference between rate of activation and inactivation with rate of activation increasing faster than the rate of channel inactivation. In Region II of Figure 16A, increasing SR Ca2+ content begins to translate into a substantial increase in peak [Ca2+] myo . This is because increasing diastolic [Ca2+] jSR causes a significant increase in the peak RyR open probability, as shown in Figure 16C, which translates into a commensurate increase in the peak of the cytosolic Ca2+ transient. The inhibiting role of the luminal sensor is relieved with building [Ca2+] jSR levels. The nonlinearity observed between the diastolic [Ca2+] jSR levels 673μ M (Figure 16B) and 734μ M (Figure 16C) at the transition between Regions I and II corresponds to the existence of a threshold characteristic for RyR release [107, 108]. As the diastolic [Ca2+] jSR levels increase beyond 734μ M, the area enclosed by the [Ca2+] ryr - O 2 ryr loop increases substantially reaching a maximum value at [Ca2+] jSR level of 948μ M not only due to a significant increase in both peak [Ca2+] ryr and peak O 2 ryr but also due to a rapid increase in rate of RyR activation resulting in a much larger difference in rate of activation and inactivation. Region III exhibits a large rapid increase in peak [Ca2+] myo due to the large RyR Ca2+ release associated with the high SR Ca2+ content. However, it is important to note that this is not a result of increasing peak RyR open probability, which begins to plateau for steady state [Ca2+] jSR values larger than 850μ M (traces 6-1 in Figure 16C). The continuing rapid increase in peak cytosolic Ca2+ concentration in region III is a combined result of: (a) increased release owing to large SR Ca2+ content; (b) large values for saturated RyR open probability supported by an increase in the area contained by each loop as seen in traces 6-1 of Figure 16C; and (c) saturated operation of the SERCA pump which acts as a predominant buffer in restoring the Ca2+ concentration levels in the cytosol after RyR release . As the diastolic [Ca2+] jSR levels increase beyond 948μ M, the area enclosed by the [Ca2+] ryr - O 2 ryr loop begins to decrease despite increasing peak [Ca2+] ryr due to saturation of the peak RyR open probability and increasing rate of RyR channel inactivation due to large values of local Ca2+ concentration ([Ca2+] ryr ) assisting in faster recovery . This model generated relationship between the peak [Ca2+] myo and the pre-release [Ca2+] jSR agrees with Figure 4 in Diaz et. al. .
As seen in traces 1-10 of Figure 16C, during the initial activation of a RyR channel (linear region of phase A in Figure 15), at any specific [Ca2+] ryr the RyR open probability takes a lower value despite a higher steady state diastolic [Ca2+] jSR level. This decreasing slope with increasing diastolic, pre-release [Ca2+] jSR reflects a faster rate of rise in [Ca2+] ryr for a larger pre-release SR content with a very small accompanying increase in the RyR open probability.
Ca2+ Release and its Effect on ICa,L
Adachi-Akahane et al.  investigated I Ca,L -induced Ca2+ release in rat ventricular myocytes in the presence and absence of caffeine, which altered [Ca2+] myo . Using their protocols and data as a guide, we modeled the effect of caffeine on jSR Ca2+ release and its subsequent effect on I Ca,L . Specifically, a Boltzman relationship was adopted to simulate the binding of caffeine to the ryanodine receptor increasing channel conductance, which in turn enhances jSR Ca2+ release. The process is assumed to be dose-dependent, with an open probability of 0.5 at a caffeine concentration of 5 mM (Appendix A1, Eq. 86).
Figures 17A and 17B show model-generated [Ca2+] myo responses in myocytes dialyzed with 2 mM Fura-2, where a simple voltage clamp pulse (-40 mV to 0 mV) is applied in the presence and absence of caffeine. I Ca,L waveforms associated with these two protocols are shown in Figure 17C. In protocol A, conducted at normal rest levels of [Ca2+] myo (10 nM), the depolarizing test pulse to 0 mV from the holding potential of -40 mV for 0.1 sec fully activates I Ca,L , which in turn triggers a rapid but small amplitude Ca2+ transient (a rise from 30 nM to 95 nM as shown in Figure 17A. Protocol B starts with an application of 5 mM caffeine for a period of 1.0 sec (open state probability of RyR channel is 0.5), followed by the same depolarizing test pulse applied at 1.0 sec. The conditioning caffeine pulse induces a [Ca2+] myo transient rising from a resting value of 30 nM to a peak value of 110 nM (Figure 17B, first peak -note the different timescale). The short depolarizing pulse activates I Ca,L , which triggers a much smaller [Ca2+] myo transient via CICR (30 nM, Figure 17B, second peak). This I Ca,L -evoked Ca2+ transient in response to the voltage pulse is smaller in amplitude than control (Figure 17A), due to: (a) reduction in RyR release due to depletion in [Ca2+] jSR ([Ca2+] jSR drops from 1.25 mM to 0.25 mM with caffeine application).; and (b) the strong competing role of the luminal RyR sensor in inactivating the RyR channel in response to the post-caffeine, low SR content. The I Ca,L -inactivation characteristics undergo a change as shown in Figure 17C, where trace B inactivates more slowly compared to the control I Ca,L (trace 1). After exposure to caffeine, the overall amplitude of [Ca2+] dhp during release does not increase as much as in the control case, since jSR Ca2+ release is reduced. Consequently, Ca2+ dependent inactivation (during release) of I Ca,L after caffeine application is not as strong as that in the control case. However, the baseline level of Ca2+ concentration at the mouth of the DHP channel ([Ca2+] dhp ) is elevated following caffeine application, causing sustained inactivation resulting in the crossover of the current trace in Figure 17C. This is consistent with the observation of Cens et al. , that Ca2+ induced inactivation is a result of calmodulin mediated sensing of the local Ca2+ concentration.
Analogous to the experimental protocol of Adachi-Akahane et al.  dealing with thapsigargin blockade of the SERCA pump (their Figure 9), we apply a series of voltage pulses (-40 mV to 0 mV, duration 0.1 sec, frequency 0.05 Hz) to our cell model, modified by blockade of the SERCA pump and addition of 2 mM Fura-2. Figure 18 shows that when Ca2+ uptake is blocked, SR Ca2+ content declines, as does the jSR Ca2+ release with each voltage pulse (Figure 18A). Considerable time is consumed after thapsigargin application, before SR Ca2+ content is finally exhausted (8-10 beats at 0.05 Hz; Figure 18A). As the magnitude of the [Ca2+] myo transient decreases, the [Ca2+] dhp dependent inactivation of I Ca,L slows (Figure 18B), resulting in an enhanced peak.
For a given clamp amplitude, a decrease in jSR Ca2+ content decreases Ca2+ release, producing decreased Ca2+ concentration at the dyadic side of the DHP-sensitive channel. Consequently, a smaller amount of calcium-calmodulin complex (Ca4CaM) is activated, leading to a lower degree of Ca2+ dependent inactivation. Thus, I Ca,L peak is enhanced, as shown in Figure 18B, increasing the area under the waveform, corresponding to an increased amount of Ca2+ entering the cell. This inverse relationship between jSR Ca2+ content and the amount of trigger Ca2+ entering the dyadic space for a given clamp voltage is the essence of the feedback mechanism maintaining the efficacy of CICR. This mechanism is operational even in the presence of high concentrations of cytoplasmic Ca2+-buffers, since the small dimensions of the dyadic space allow jSR released Ca2+ to reach its destination (i.e., the DHP receptor) before being captured by the buffers .
Effect of modulation of [Ca2+]o
Increase in [Ca2+] o causes an increase in entry of extracellular Ca2+ into the cell via the I Ca,L channel. This results in an increase in SR Ca2+ content, as shown in Figure 19B. Aided by an increase in trigger Ca2+ as well as a larger SR Ca2+ content, the RyR release is enhanced. This increase in SR release assisted by impaired Na+/Ca2+ exchange due to elevated [Ca2+] o manifests in an increase in the peak Ca2+ concentration levels in the cytosol, as shown in Figure 19A. Similarly, a decrease in extracellular Ca2+ levels causes a decrease in trigger Ca2+ available to cause CICR, which results in a decrease in peak [Ca2+] levels in the cytosol, as shown in Figure 19C. Reduced Ca2+ entry into the cell via I Ca,L , translating into decreased availability of post-release Ca2+ in the cytosol, causes reduction in SR filling via the bi-directional SERCA pump. The corresponding drop in the SR Ca2+ content is evident in Figure 19D. Although, our model-generated results show similarity to data in Figure 2, Diaz et al. , our result in Figure 19D does not confirm an increase in SR Ca2+ content with a decrease in [Ca2+] o as reported in Figure 2B, Diaz et al. . On the contrary, they agree with the data in Figure 5, Diaz et al. 1997 which shows a decrease in free diastolic steady-state SR Ca2+ content with decrease in [Ca2+] o in cells not showing spontaneous release. In this study, we have attempted to isolate the effects of changes in [Ca2+] o on intracellular Ca2+ levels by ensuring that there are no voltage dependent effects on Ca2+ transport via the I Ca,L channel, Na+/Ca2+ exchanger or the plasma membrane Ca2+ ATPase pump. The simulation protocol used in our study (reproduced from Diaz et al. ),is a pulse train of amplitude (-40 mv to 0 mv), duration (100 ms) and a frequency of 0.5 Hz. However, it is important to note that, for any constant extracellular Ca2+ concentration, in response to a voltage clamp pulse train, increasing clamp potential (> 10 mv) results in a decrease in steady state SR Ca2+ content as a result of the dominating effect of reduced Ca2+ entry via the I Ca,L channel despite a decrease in net Ca2+ extrusion via the Na+/Ca2+ exchanger per cycle.
Calcium balance under conditions of repetitive stimulation
Our model exhibits long term Ca2+ stability. Analogous to the experimental protocols of Negretti el al. , we studied the dynamic aspects of the calcium balance in the cell model by subjecting it to repetitive voltage clamps of different durations. For stable steady-state operation, Ca2+ entry into the cytosol via I Ca,L must exactly balance Ca2+ efflux. Changing the rate or pattern of stimulation can have significant effects on the cell's cytosolic Ca2+ balance and subsequent contractile response . To demonstrate Ca2+ balance in our model, long (Figure 20A) or short (Figure 20B, C and 20D) voltage clamp pulses were applied at a selected repetition frequency, and the resultant changes in [Ca2+] myo were correlated with sarcolemmal and SR Ca2+ fluxes. Analogous to the protocol of Negretti et al. , pulse trains were applied from a holding potential of -40 mV to +10 mV, with durations of either 800 ms (long-pulse) or 100 ms (short-pulse) at a frequency of 0.33 Hz, for a total test duration of 2 minutes (40 pulses). In analyzing the Ca2+ balance, the amount of electric charge Q transferred per pulse was calculated according to the equation: , where I(t) can be I Ca,L , I NaCa , I PMCA , I ryr , I cyt,serca , and T is 1 cycle duration (3 ms for 0.33 Hz stimulation).
Numerical integration was carried out with respect to steady-state levels, hence only activating currents and leakage currents were considered (background Ca2+ current is neglected). The results of the Ca2+ balance are presented in a manner where relative Ca2+ fluxes can be easily compared. We note that in evaluating the integral of the combined exchanger current I NaCa , the integral is multiplied by a factor of 2 to account for the fact that the exchanger stoichiometry is 3Na+ per Ca2+ (the exchanger transports 1 net charge per Ca2+, whereas all other Ca2+ currents transport 2 charges per Ca2+).
During repetitive long-pulse (800 ms) stimulation, transient Ca2+ fluxes cross the sarcolemmal and SR membranes in both directions (inward (into the cytoplasm) and outward). However, [Ca2+] myo returns to control levels by the end of the long-pulse stimulation protocol. The long-pulse stimulation protocol thus serves as a means of studying the steady-state balance of Ca2+ influx and efflux to and from the cytoplasm, under control conditions. In contrast, this Ca2+ balance is not present in the short-pulse protocol, and [Ca2+] myo is elevated at pulse termination. Our long-pulse simulations show that: (1) the magnitude of individual Ca2+ transients do not change during the 40 pulse sequence; (2) the peak calcium release current (I ryr ) has a constant magnitude; and (3) the occupancy of the calsequestrin Ca2+ buffer in the jSR compartment is reduced by 15% as a result of Ca2+ release during each cycle. Although Ca2+ fluxes cross the sarcolemma and SR membranes in either direction, the integral of each of these currents (indicating charge transfer) over the pulse repetition interval is a constant. Figure 20A shows this well, in that, the charge transfer for each model current (large or small) is a constant, indicating no net loss of Ca2+ nor consecutive-pulse Ca2+ depletion (or augmentation) in the jSR compartment under the long-pulse protocol. In Figure 20, we refer to Ca2+ influx and efflux through the bounding cell membrane (sarcolemma (SL) and T-tubule) as "external fluxes" (i.e., I Ca,L , I NaCa , and I PMCA ), whereas SR membrane Ca2+ fluxes (I ryr , I cyt,serca ) are called "internal fluxes." By convention, inward fluxes are negative (e.g. I Ca,L ) and outward positive (e.g., I NaCa and I PMCA ). Figure 20A shows that the integrated external fluxes sum to zero. The average Ca2+ charge entering the cell via I Ca,L is 603.28 pC (pico-coulomb), and the sum of averaged Ca2+ charges extruded from the cell via I NaCa and I PMCA is 603.28 pC (602.2 pC by I NaCa and 1.08 pC by I PMCA . Figure 20A also shows the magnitudes of the integrated internal fluxes, which sum to zero as well. In calculating inward and outward Ca2+ fluxes relative to the SR lumen, we consider an inward flux (e.g., I cyt,serca ) negative and an outward flux positive (e.g., I ryr ). The integral of I ryr is 751.8 pC, compared with integral of the I Ca,L trigger current (603.28 pC). The calcium gain of CICR calculated in this way is low owing to the increased entry of Ca2+ via the trigger current I Ca,L for a prolonged duration of 800 ms compared to 50 ms previously. The integral of SR Ca2+ uptake is 751.8 pC. These two internal component currents yield flat integral values during long-pulse stimulations and sum to zero.
Decreasing the length of the depolarizing pulse to 100 ms (short pulse) has a pronounced effect on the recovery of contraction (Figure 20B, C and 20D). Compared with the relatively flat amplitude of [Ca2+] myo observed with stimulation by long pulses, short pulse stimulation produces a negative staircase (Figure 20B) in agreement with the data (Negretti et al. ), wherein the peak magnitude decays exponentially to a steady state level of 60% of the initial value. Figure 20B also shows that the model predicts a similar exponential decay in jSR Ca2+ concentration [Ca2+] jSR . Both of these indicators show a net loss of Ca2+ from the cell. Figure 20C and 20D show the transient charge movements associated with the external and internal Ca2+ currents respectively during application of the short-pulse train, which are indicative of an elevated Ca2+ load and resultant Ca2+ imbalance when the pulse train is first applied. Balance is achieved, but at new lower values of [Ca2+] myo and [Ca2+] jSR . For convenient comparison, all internal fluxes are made positive and all external fluxes negative. Figure 20C shows that I NaCa plays a significant role in forming the transient response seen in the short-pulse protocol in extrusion of excess Ca2+ load on the cell. The decreasing amplitude of the Ca2+ transient causes a decreased LSR uptake current (I cyt,serca ), which in turn causes a decrease in [Ca2+] jSR resulting in decreased release (I ryr ). With successive pulses, I cyt,serca and I ryr decline further, mirroring the exponential decline in I NaCa toward balance. At steady state, the external and internal fluxes reach the new equilibrium, where the sum of the integrals of external currents and the sum of the integrals of internal currents are zero.
Long-term calcium stability at higher pacing rates
Most VC experiments on cardiac cells are usually conducted at low stimulation rates. Our simulations of typical VC experiments at these low rates all exhibit long term calcium stability (Figures 17, 18, 19 and 20). However, our model can also exhibit sustained calcium balance at higher pacing rates. Figure 14A shows steady state cytosolic Ca2+ transients in response to a repetitive 4 Hz (which is more physiological for a rat ventricular myocyte) voltage clamp stimulation lasting for 2.5 seconds. The sustained peak and basal Ca2+ levels of [Ca2+] myo over this prolonged time frame are indicative of long tern Ca2+ stability in the model. The corresponding [Ca2+] jSR profile over 2.5 seconds is also shown in Figure 14A, which indicates sustained SR filling owing to the luminal sensor mediated control of the RyR channel. The stimulation protocol used is a pulse train of amplitude (-40 mv to 10 mv), duration (50 ms) and frequency of 4.0 Hz. Tables 9 and 10 provide values for the initial conditions used in our model.
Secondary [Ca2+]myo transients induced by "tail currents"
Our voltage clamp simulations have dealt with clamp voltages in the range -30 mV ≤ V ≤ 40 mV. Some experimental studies [100, 2, 112] however have employed even larger clamp voltages (40 mV ≤ V ≤ 60 mV) to explore the high voltage behavior of the DHP-sensitive calcium channel. Such large clamp voltages elicit a brief Ca2+ influx called a "tail current", which has been shown to trigger RyR release and hence cause contraction during repolarization [100, 2]. The secondary Ca2+ transient induced by the "tail current" is a critical argument in favor of CICR as these "tail currents" are not observed in skeletal muscle where the membrane potential directly controls SR Ca2+ release . Our model reproduces this secondary Ca2+ transient observed during the return to resting potential from a large clamp voltage (≥ 40 mV). Figure 21A shows a cartoon depicting the voltage clamp stimulation protocol used in our study where a pulse of amplitude (-40 mV to +50 mV) is employed with the pulse duration (T p ; Figure 21A) increasing from 50 ms (trace 1) to 200 ms (trace 7) in steps of 25 ms. The clamp potential transition time (T t ) is fixed at 1 ms. The peak of the "tail current" at the end of the pulse is 25 times larger than the peak of the I Ca,L current observed at the beginning of the pulse. This is in agreement with model generated VC data reported by Geenstein et al.(; Figure 8). Figure 21B shows that, with an increase in clamp pulse duration, there is a corresponding increase in the peak of the trail current (see traces 1-7) until the effect ultimately saturates. Corresponding to this increase in peak I Ca,L tail current, there is an increase in the open probability of the RyR Ca2+ channel, with RyR Ca2+ release indicated by the corresponding decrease in [Ca2+] jSR (panel C) and increase in [Ca2+] myo (panel D). We note from Figure 21C that for a shorter 50 ms VC pulse, there is no RyR Ca2+ release although an I Ca,L tail current is produced. This tail current produced with the 50 ms pulse slightly augments the [Ca2+] myo transient just after its peak (trace 1, panel E), but there is no corresponding decrease in the [Ca2+] jSR transient (trace 1, panel D) indicating no RyR Ca2+ release. At the end of this short pulse (50 ms), the RyR channel is in fact in its absolute refractory period. The small secondary increase in [Ca2+] myo transient as seen in data corresponding to trace 4 of Figure 7C is a result of the large Ca2+ influx via the I Ca,L "tail current". During this phase, the Na+/Ca2+ exchanger is biased to extrude Ca2+ out of the dyad, and hence cannot be responsible for this secondary rise in the cytosolic Ca2+ transient. As pulse duration is increased beyond 100 ms the "tail current" causes a gradual increase in RyR open probability (traces 3-7 in Figure 21) indicating progressive recovery from the refractory period. Minimum value of [Ca2+] jSR in traces 3-7 of Figure 21D show increase in SR Ca2+ release with increasing pulse duration as a result of adequate RyR channel recovery and increasing I Ca,L "tail current". Cytosolic Ca2+ transients shown in Figure 21E indicate "tail current" induced secondary RyR release for pulse duration ≥ 100 ms. It has also been previously reported  that recruitment of an additional population of previously 'silent' Ca2+ channels could cause facilitation of tail currents at increasingly large clamp voltages with a time-dependence associated with the recruitment process. Our model not only accounts for the contribution of the already open channels (in states O 2 dhpr & O 3 dhpr ) to the tail current during repolarization but also allows for the voltage- and time-dependent recruitment of 'silent' Ca2+ channels modeled using a high voltage state C 6 dhpr . This study proves to be an adjunct to the study presented in Figure 13 in establishing the refractory nature of the RyR channel.
In our model, cytosolic buffering is attributed to several factors including: (a) calmodulin (CaM); (b) the Ca2+-specific (Tc) troponin binding site; (c) the Ca2+ - Mg2+ competitive troponin binding site; and (d) the fluorescent indicator dye Fluo3 used to detect changes in [Ca2+] myo . The effects of the component buffers in helping to maintain a low [Ca2+] myo are shown in Figure 22 in terms of occupancy functions, such as O C (fractional occupancy of calmodulin by Ca2+), O tc (fractional occupancy of troponin-Ca sites by Ca2+), O tmgc (fractional ocupancy of troponin-Mg sites by Ca2+), O tmgmg (fractional occupancy of troponin-Mg sites by Mg2+) in the cytosol, and O calse (fractional occupancy of calsequestrin by Ca2+) in the jSR Ca2+ release compartment. Figure 22 shows that when the Ry-sensitive Ca2+ release channel is triggered, the jSR releases Ca2+ and the occupancy O calse in the jSR, declines from 68% to 37%. Ca2+ release from the jSR induces fast Ca2+ binding by calmodulin and troponin in the cytosol, as represented by the increases of O c (0.29) and O tc (0.16) in Figure 22. Interactions between Ca2+ and the troponin-Mg sites result in an increase in occupancy of these sites by Ca2+ (O tmgc ), and a decrease in occupancy by Mg2+ (O tmgmg ).
It is well-established that mammalian cardiac excitation-contraction coupling is mediated by calcium-induced calcium release (CICR). We have developed a comprehensive mechanistic model of CICR under voltage clamp conditions in the rat ventricular myocyte, which includes electrical equivalent circuit models for both the free sarcolemma and that portion involving junctional transmission, as well as, fluid compartment models for several fluid media within the cell (dyadic cleft space, longitudinal sarcoplasmic reticulum (LSR or Ca2+ uptake compartment), junctional sarcoplasmic reticulum (jSR or release compartment), and the cytosolic fluid compartment). An external bathing medium completes our fluid compartment description of the cell Figure 1. The multiple component model is referred to as the "whole-cell model" (Figures 1 and 2). We have probed the mechanisms regulating Ca2+ in the cell. In particular, we have focused on the dyadic mechanisms effecting CICR.
The dyadic controller is a finely tuned coupling device consisting of two opposed Ca2+ channels separated by a small dyadic space (Figure 1B): the sarcolemmal DHP-sensitive "trigger" channel and the Ry-sensitive jSR Ca2+ release channel. The trigger channel is voltage-activated and is driven in our simulations by a voltage clamp pulse which opens the channels, admitting Ca2+ influx (trigger current) to the dyadic space. After diffusion within the dyadic space, trigger Ca2+ affects the Ca2+ dependent open probability of the jSR Ca2+ release channel over a small range of Ca2+ concentrations. Clamp voltage magnitude strongly affects the I Ca,L current and hence the amount of Ca2+ delivered to the dyadic space, but the range over which the Ca2+ concentration at the mouth of the Ca2+ release channel changes (in the absence of RyR release) is quite small (23.42μ M to 8.82μ M for a change in VC from 10 mv to 40 mv). The RyR-mediated Ca2+ release is activated by trigger current, but the release itself is affected by the concentration gradient between the jSR and the dyadic space, and the temporal and refractory properties of the ryanodine receptor (Figure 13). The four-state RyR model is informed regarding supply Ca2+ content by the jSR luminal sensor, a novel feature in our model which characterizes the important protein-protein interactions between calsequesterin, triadin/junctin and the RyR receptor. Triadin/junctin strongly regulates the sensitivity of RyR to trigger Ca2+. Hence the luminal sensor, which is a key element responsible for robust post-release RyR inactivation and refractoriness of the Ry-sensitive Ca2+ release channel, is critical in providing realistic fits to cytosolic Ca2+ transients and an adequate refilling time for the SR Ca2+ stores. Ca2+ release is thus graded with Ca2+ concentration at the mouth of Ca2+ channel in a very sensitive manner with a gain of approximately 7, as is shown in Figure 12.
The sarcolemmal portion of the dyadic membrane defining the dyadic space contains voltage-sensitive Ca2+ channels and deals with changes in the external environment of the ventricular cell (e.g., membrane response to changes in transmembrane potential and chemical signaling agents). It integrates these various stimuli and delivers a trigger current to the small dyadic space. In contrast, the jSR membrane lining the opposite boundary of the dyadic space is concerned with the adequacy of Ca2+ release in CICR. It contains Ry-sensitive Ca2+ channels that require a Ca2+ concentration gradient directed across the channel and into the dyadic space for their operation. Thus, jSR Ca2+ concentration must be maintained within an acceptable range so that calcium is always available for ready release. The relationship between jSR Ca2+ content and peak [Ca2+] myo is shown in Figure 16A. The LSR compartment connects the jSR compartment with the cytosolic compartment. It feeds make-up calcium to the jSR, by using a SERCA pump to actively draw in Ca2+ from the myoplasm. Pumping rate is controlled by the cytosolic Ca2+ concentration [Ca2+] myo as well as the LSR Ca2+ concentration [Ca2+] LSR as part of a mechanism for replenishing and maintaining jSR Ca2+ stores. Our model also incorporates the Ca2+ induced CaM mediated effects of CaMKII and CaN on targets such as the DHP-sensitive I Ca,L channel, Ry-sensitive Ca2+ release channel, as well as the SR Ca2+ ATPase pump. Provision for the effects of phospholamban on SERCA has also been included.
The voltage-gated and chemically gated channels of the dyad are tightly coupled by feedback mechanisms that involve Ca2+ signaling. Although physically separated, the voltage-sensitive Ca2+ channel (Figure 9) as well as the Ry-sensitive release channel (Phase C in Figure 15) are inhibited by increased Ca2+ levels in the dyadic space . In earlier models which lacked the luminal sensor, RyR self-inhibition by dyadic Ca2+ was the only mechanism besides stochastic attrition (both of which are inadequate in effecting a robust RyR channel closure) that was given the role of initiating RyR recovery. This has been tested in our model by artificially clamping all the variables (including concentration at the mouth of the RyR channel) to the levels reached at T2 in Figure 15 despite which the RyR open probability, begins to decline in phase C of Figure 15 showing the self-inhibitory role of high dyadic Ca2+ concentration. Our model predicts (data not shown) that this Ca2+ signaling continues even in the presence of high concentrations of Ca2+ buffering agents in the cytosol (in agreement with the data of Adachi-Akahane et al.  and Diaz et al. ). Tight coupling between the DHP and Ry-sensitive Ca2+ channels within the dyadic space thus preserves the mechanism of CICR under these extreme conditions.
Regulation of cytosolic Ca2+ concentration [Ca2+] myo is evident at the cellular level. It is affected strongly by the voltage and [Ca2+] myo -dependent properties of the sarcolemmal currents I NaCa and I Ca,L , as well as the [Ca2+] myo -dependent properties of the I PMCA and SERCA pumps. We demonstrate a model-generated whole-cell Ca2+ balance, which shows the importance of the Na+/Ca+ exchanger in extruding the Ca2+ that has entered the cell under normal activity, and also any excess that might occur when cytosolic Ca2+ levels rise. The variety of experiments emulated in this study demonstrates quantification of Ca2+ balances for all external and internal Ca2+ fluxes and shows that the model has long-term stability in regulating cytosolic Ca2+, as shown in the 120-sec duration experiments of Negretti et al. (Figure 20) at a pulse repetition rate of 0.33 Hz., and the faster paced stimulation at 4 Hz shown in Figure 14A.
We have examined the dyadic component of the I Ca,L with regard to its Ca2+ -dependent inactivation as a function of jSR Ca2+ content (Figure 17C). Trigger current is voltage-dependent and therefore parameterized by clamp pulse amplitude (Figure 7A). In addition, its inactivation is Ca2+ dependent, especially during Ca2+ release. At a constant depolarizing voltage pulse, the pre-release jSR Ca2+ content dictates the peak of the cytosolic Ca2+ transient (Figure 16A), by controlling peak Ca2+ concentration in individual dyads. As jSR Ca2+ content increases, so does the peak dyadic Ca2+ concentration, and the amount of sarcolemmal Ca2+ influx declines due to greater Ca2+ dependent inactivation. This autoregulatory feedback mechanism helps to establish a stable operating point for jSR Ca2+ content. Transient increases in jSR Ca2+ content bring about increases in Ca2+ release, but reflexly, decrease sarcolemmal Ca2+ influx via I Ca,L . Coupled with other dyadic and extra-dyadic mechanisms (I NaCa ) that decrease [Ca2+] myo and hence Ca2+-uptake via the SERCA pump, jSR Ca2+ content decreases. For decreases in jSR Ca2+ content, the opposite occurs, so that regardless of the sign of the perturbation in jSR Ca2+ content, it tends to stay constant in the steady-state. This is an important feature of the dyadic mechanism that preserves the integrity of CICR. The gain of this feedback system about the operating point is visualized as the slope of the peak Ca2+ transient vs SR Ca2+ content characteristic (Figure 16A) for a given voltage pulse level. Observing this figure, we note that there is a linear operating range beyond which the system gain increases dramatically. Our model results also show that a decrease of jSR Ca2+ content (simulated either by a decrease of Ca2+ uptake into the LSR by thapsigargin or an increase of Ca2+ leak out of the jSR by caffeine) decreases systolic [Ca2+] myo , and hence the model might serve as a useful adjunct in a study of heart failure, where decreased contractility as a result of diminished Ca2+ transients are commonly observed .
This model of a rat ventricular myocyte is limited to Ca2+ related channel, exchanger and pumps (I Ca,L , I NaCa , I PMCA and SERCA pump), while lacking exclusive Na+ or K+ related channels and transporters and is based on data at positive potentials in the range 10 mV ≤ V ≤ 40 mV. It is aimed at mimicking voltage clamp conditions where channels other than calcium are blocked, and it cannot be used to study experiments involving the generation of action potentials. However, its focus on the Ca2+ dynamics allows one to comprehend more clearly the important role of Ca2+ signalling pathways and feedback control systems in maintaining whole cell homeostasis over a prolonged period of time.
A single dyadic space in our model has one representative, lumped DHP-sensitive and Ry-sensitive Ca2+ channel, on opposing sarcolemmal and SR surfaces respectively. This simplified configuration does not conform to detailed structural information regarding the geometrical relationships between DHP and RyR-sensitive Ca2+ channels [115, 55, 116] and thus cannot be used to draw conclusions about this part of the EC coupling process. However, our lumped abstraction forms a functional model of the dyadic coupling unit that can produce accurate predictions of cytosolic Ca2+ transients. The effectiveness of this model is further demonstrated by its ability to accurately characterize the interaction between the DHP and Ry-sensitive Ca2+ channels, including pulse duration dependent termination of release (Figure 11), Ca2+ dependent inactivation of I Ca,L (Figures 7, 9), as well as the wide variety of whole-cell voltage clamp protocols (Figures 17, 18, 19 and 20).
Although our model provides secondary Ca2+ tail transients elicited by I Ca,L "tail currents" (Figure 21), this aspect of the model has not been verified extensively due to paucity of measured VC data showing tail transients over the high voltage range (40 ≤ V ≤ 60) in rat ventricular myocytes. Our I Ca,L tail currents are in general agreement with model generated data at a clamp voltage of 50 mV reported by Greenstein et al. (; Figure 8). The restitution time for the RyR channel in rat ventricular myocytes is believed to be at least 25 ms  and as fast as 150 ms  indicating the wide range of values reported in the literature. We demonstrate model-generated RyR refractoriness by changing the duration of the simulated VC pulse and obtain RyR recovery characteristics that are consistent with measured data for Ca2+ spark restitution in rat ventricular myocytes reported by Sobie et al. . However, our results on this important phenomena, particularly the onset of tail transients with increasing inter-stimuli interval are preliminary and further modeling investigations would benefit considerably from the availability of additional measured data.
We have developed a mathematical model of Ca2+ dynamics under voltage clamp conditions in the rat ventricular myocyte, which is based solidly on experimental data and includes the most extensive description available of a novel feature, namely the luminal Ca2+ sensor in the junctional SR which models the protein-protein interaction between triadin/junctin, calsequestrin and the RyR channel. The luminal sensor imparts the much needed refractoriness to the Ry-sensitive Ca2+ release channel. This element is critical in providing realistic fits to cytosolic Ca2+ transients and an adequate refilling time for the SR Ca2+ stores. Our voltage-clamp simulations demonstrate graded Ca2+ transients with sufficient gain, as well as quantification of Ca2+ balances for all external and internal Ca2+ fluxes. Our model of the dyadic coupling unit (DCU) provides mechanistic explanations of the major input-output relationship for CICR (Figure 16), as well as its modulation by trigger current (clamp voltage). The variety of experiments emulated in this study demonstrates that the model has long-term stability in regulating cytosolic Ca2+, as shown in the 120-sec duration experiments of Negretti et al. (Figure 20) at a pulse repetition rate of 0.33 Hz., and the faster (physiological) paced stimulation at 4 Hz shown in Figure 14A. It also provides biophysically based insights into the molecular mechanisms underlying whole-cell responses to the wide variety of testing approaches used in voltage clamp studies of myocytes that have appeared in the literature over the past two decades (Figures 17, 18, 19 and 20). Thus, the model serves as a platform for the predictive modeling of VC investigations in a number of areas. These include new hypotheses with regards to the under-expression of triadin/junctin resulting in a malfunctioning luminal sensor, which could affect long-term calcium stability of the cell (Figure 14B), and/or changes in the refractoriness of the RyR Ca2+ channel (Figures 13 and 21) affecting the integrity of CICR under a variety of conditions. These are fundamental issues that would benefit from a better mechanistic understanding of deranged calcium signalling in the rat ventricluar myocyte. This study is aimed at providing an initial step towards this goal.
Below is the complete set of equations used in the model.
A1 - Equations for currents in the model
L-Type Ca 2+ current
Ca2+ current through the DHP-sensitive ICa,Lchannel
Na+ current through the DHP-sensitive ICa,Lchannel
Cs+ current through the DHP-sensitive ICa,Lchannel
where k Cs = 0.5 and [Ca2+] dhp , [Na+] dhp , [Cs+] dhp are Concentrations at the mouth of the DHP-sensitive Ca2+ channel.
The corresponding unitary currents i Ca , L , i Na , i Cs are obtained by dividing the above net channel currents by the number of dyadic units, N dyad .
Gating scheme for the 2-state Markovian model used to allow Ca2+ mediated interaction of the DHP-sensitive I Ca , L channel and calmodilin:
Expressions for the rate constants:
Gating scheme for the 6-state Markovian model for the DHP-sensitive L-type Ca2+ release channel:
The open state O 3 dhpr accounts for the increased tail current produced as the result of a large depolarization. Expressions for rate constants:
• Case 1: With CICR (control)
•Case 2: With Ryanodine applied (no RyR release)
•Case 3: With Ca2+ substituted with Ba2+
Uptake of Ca2+ from the cytosol into the LSR
Ca2+ fluxes from cytosol to SERCA and SERCA to LSR
Differential equation for Ca2+ buffered by the SERCA protein:
Expressions for the rate constants for Ca2+ binding to/release from SERCA:
Affinities for the forward and backward Ca2+ fluxes:
Differential equation for phospholamban
Ca2+ pump in SL
where [Ca2+] NaCa , [Na+] NaCa are Concentrations at the mouth of the NaCa-exchanger.
The corresponding unitary current i NaCa is obtained by dividing the above net channel current I NaCa by the number of dyadic units N dyad .
The corresponding unitary current i NaCs is obtained by dividing the above net channel current I NaCs by the number of dyadic units N dyad .
Background Na+ current
Ca2+ transfer from LSR to a single jSR
Ca2+ release from a unit jSR into a single DCU
Gating scheme for the 4-state Markovian model for the RyR-sensitive SR Ca2+ release channel:
where [Ca2+] ryr is the Ca2+ concentration at the mouth of the RyR channel on the dyadic side. In the presence of caffeine (CF, concentration in μ M):
Gating scheme for the 6-state Markovian model for the Luminal Calcium sensor.
A 1 ls Fractional occupancy of RyR by Triadin/Junctin
I 2 ls Fractional occupancy of RyR by Triadin/Junctin and Calsequestrin
I 3 ls Fractional occupancy of Triadin/Junctin by Calsequestrin
I 4 ls Free Triadin/Junctin, B5ls - Free Calsequestrin
B 6 ls Fractional Occupancy of Calsequestrin by Calcium
Ca2+ dependent CaM mediated activation of CaMKII and CaN
Ca2+ binding to CaM and CaM buffering
A2 - Differential equations for buffers used in the model
Fluorescent indicator dye
Intracellular Ca2+ buffering:
(Fractional occupancy of troponin-Ca complex by Ca2+):
(Fractional occupancy of troponin-Mg complex by Ca2+):
(Fractional occupancy of troponin-Mg complex by Mg2+):
A3 - Differential equations for ion concentrations used in the model
Intracellular Ca 2+ concentration:
1. Ca2+ concentration in the cytosol
where I dyad is the net integrated Ca2+ flux diffusing out of all the dyadic units into the cytosol
2. Ca2+ concentration in the jSR
3. Ca 2+ concentration in the LSR
4. Diffusion equation for Calcium in the dyadic space
Intracellular Na + concentration:
1.Na+ concentration in the cytosol
2. Na+ concentration in the dyadic space
Intracellular Cs + concentration:
1. Cs+ concentration in the cytosol
2. Cs+ concentration in the dyadic space
calcium ion concentration
- [Ca2+] dhp Ca:
2+ concentration at the mouth of the DHP-sensitive I Ca,L channel
- [Ca2+] dyad spatial Ca:
2+ concentration in the dyad
- [Ca2+] jSR luminal Ca:
2+ concentration in the jSR
- [Ca2+] LSR Ca:
2+ concentration in the LSR
- [Ca2+] myo cytosolic Ca:
- [Ca2+] serca Ca:
2+ concentration buffered by the SERCA protein
- [Ca2+] o extracellular Ca:
- [Ca2+] ryr Ca:
2+ concentration at the "mouth" of the RyR channel on the dyadic side
Ca2+ bound calmodulin
Ca2+ bound to the buffer Fluo3
Ca2+/calmodulin-dependent protein kinase II
- CaMKII act activated Ca:
2+/calmodulin-dependent protein kinase II
- C a N act :
calcium dependent facilitation
calcium dependent inactivation
- [Cs+] myo Cs:
+ concentration in the cytosol
- [Cs+] dyad Cs:
+ concentration in the dyadic space
dyadic coupling unit
affinity of backward Ca2+ flux from LSR to cytosol
affinity of forward Ca2+ flux from cytosol to LSR
- I NaCa :
sodium calcium exchanger current
- I NaCs :
sodium cesium pump current
- I PMCA plasma membrane Ca:
2+ ATPase pump current
- I ryr Ca :
2+ current due to CICR from an individual jSR
- I tr Ca :
2+ current due to concentration gradient driven Ca2+ transport from LSR to jSR
junctional portion of the sarcoplasmic reticulum
- k d :
- k mp half saturation constant for the sarcolemmal Ca:
L-type DHP-sensitive Ca2+ channel
long lasting type
longitudinal portion of the sarcoplasmic reticulum
- [Na+] dyad Na:
+ concentration in the dyadic space
- [Na+] myo Na:
+ concentration in the cytosol
- O c fractional occupancy of calmodulin by Ca:
2+ in the cytosol
- O calse fractional occupancy of calsequestrin by Ca:
2+ in the jSR
- O tc fractional occupancy of troponin-Ca sites by Ca:
2+ in the cytosol
- O tmgc fractional ocupancy of troponin-Mg sites by Ca:
2+ in the cytosol
- O tmgmg fractional occupancy of troponin-Mg sites by Mg:
2+ in the cytosol
- P o :
protein kinase A
- PLB dp :
- PLB p :
phospholamban to SERCA ratio
sarcoplasmic reticulum Ca2+ ATPase
Ca2+-specific troponin binding site
voltage dependent inactivation
Fabiato A, Fabiato F: Contractions induced by a calcium-triggered release of calcium from the sarcoplasmic reticulum of single skinned cardiac cells. J Physiol. 1975, 249: 469-495.
Beuckelmann DJ, Wier WG: Mechanism of release of calcium from sarcoplasmic reticulum of guinea-pig cardiac cells. J Physiol(Lond). 1988, 405: 233-235.
Stern MD: Theory of excitation-contraction coupling in cardiac muscle. Biophys J. 1992, 63: 497-517. 10.1016/S0006-3495(92)81615-6.
Hilgemann DW, Noble D: Excitation-contraction coupling and extracellular calcium transients in rabbit atrium: reconstruction of basic cellular mechanisms. Proc R Soc Lond B. 1987, 230: 163-205. 10.1098/rspb.1987.0015.
Isenberg G: Cardiac excitation-contraction coupling: from global to microscopic models. Physiology and Pathophysiology of the Heart. Edited by: Sperelakis N. 1995, Kluwer Academic, 289-307.
Stern MD, Song LS, Cheng H, Sham JSK, Yang HT, Boheler KR, Rios E: Local control models of cardiac excitation-contraction coupling. J Gen Physiol. 1999, 113: 469-489. 10.1085/jgp.113.3.469.
Rice JJ, Jafri S, Winslow RL: Modeling gain and gradedness of Ca2+ release in the functional unit of the cardiac diadic space. Biophys J. 1999, 77: 1871-1884. 10.1016/S0006-3495(99)77030-X.
Greenstein JL, Winslow RL: An integrative model of the cardiac ventricular myocyte incorporating local control of Ca2+ release. Biophys J. 2002, 83: 2918-2945. 10.1016/S0006-3495(02)75301-0.