Modeling CICR in rat ventricular myocytes: voltage clamp studies

Background 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). Methods 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. Results 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. Conclusions 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.


Background
Contraction of cardiac muscle is triggered by a transient rise in intracellular Ca 2+ concentration [Ca 2+ ] myo . Sarcolemmal (SL) membrane depolarization triggers Ca 2+ influx from the extracellular medium by opening dihydropyridine (DHP)-sensitive L-type Ca 2+ channels. Following diffusion across a small sub-membrane dyadic space, this influx activates ryanodine receptors (RyRs) controlling ryanodine-sensitive Ca 2+ release channels in the junctional portion of the sarcoplasmic reticulum (jSR). Fabiato and Fabiato [1] named the process calcium-induced calcium release (CICR). Ca 2+ subsequently diffuses from the dyadic space into the cytosol. Ultimately, intracellular Ca 2+ concentration [Ca 2+ ] myo is returned to resting levels by combination of: (a) Ca 2+ buffering in the dyadic space and cytosol; (b) sequestration of Ca 2+ by sarcoplasmic/endoplasmic reticulum Ca 2+ -ATPase (SERCA)-type calcium pumps lining the longitudinal portion of the sarcoplasmic reticulum (LSR); and (c) Ca 2+ extrusion from the cytosol by Na + /Ca 2+ exchangers and Ca 2+ -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 Ca 2+ release is proportional to the influx of trigger Ca 2+ [2], whereas high gain indicates that the SL trigger current elicits a high SR Ca 2+ release flux. Graded Ca 2+ release with high gain is somewhat paradoxical according to Stern [3], 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 Ca 2+ current. Stern [3] proposed that such a gradation paradox might be explained if the stimulus for Ca 2+ release by RyRs were actually the local nanodomains of [Ca 2+ ] generated by nearby L-type channels, rather than the global cytosolic [Ca 2+ ] myo . According to this hypothesis, graded control of macroscopic SR Ca 2+ release can be achieved by graded statistical recruitment of individual, autonomous, all-or-none stochastic release events [6]. 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 Ca 2+ release complex, including several [7][8][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. [10] 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 Ca 2+ entry via the I Ca,L channel to triggered SR Ca 2+ release is known to increase with decrease in the magnitude of I Ca,L,TT [13], the modeling of which is made possible by considering Ca 2+ 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 [14] points towards the important role of the Ca 2+ refilling rate from the network to junctional SR in controlling RyR release termination via the luminal sensor. Shiferaw et al. [15] developed a computationally tractable model of Ca 2+ cycling to represent the release of calcium from the SR as a sum of spatially localized events that correspond to Ca 2+ sparks, assuming the recruitment rate of Ca 2+ 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 Ca 2+ release complex that includes: (a) a sub-sarcolemmal dyadic cleft space separating the SL and jSR membranes; (b) a single DHP-sensitive Ca 2+ 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 [16], we further assume that each ventricular cell contains 10,000 of these dyadic Ca 2+ 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 [3] local domain concept by considering the aforementioned local nanodomains identical, but focusing on the nonlinear dynamics of the two different types of Ca 2+ 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.

Experimental Methods
Rat ventricular myocytes were prepared from 200-300 g male Sprague Dawley rats by dissociation with collagenase, as previously described [17]. 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 [17]. 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 Ca 2+ current rundown. External solution in the bath was normal Tyrode (1 mM Ca 2+ ) 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.

Computational Aspects
All simulations and analysis were performed on a 2.8 GHz Intel® Core™2 Duo CPUbased 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 leastsquares method [18] was used to obtain the solution of the system of non-linear ordinary differential equations. Specifically, we have employed an algorithm given by Lau [19]. 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 Ca 2+ diffusion in the cleft space. Detailed numerical methods are similar to those presented by Smith et al. [22]. The results were visualized using Matlab by Mathworks and Origin by Microcal Software.

Model Development
Our objective was to develop a model of the rat ventricular cell which could be used to explain Ca 2+ signaling at the nanoscale level of the dyad and integrate the contributions of many dyads to produce a Ca 2+ transient and continuous Ca 2+ balance at the whole-cell level. Therefore, we start with a broad discussion of the elements of the DCU and its Ca 2+ supply (the jSR), and continue with a progressively more detailed description of the whole cell model. It is important to note that all Ca 2+ concentrations discussed in the model pertain to unbound Ca 2+ unless specified.

Membrane Classification
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 [16], 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 [16].

Channel and Exchanger Distribution
Recent research has also shown that besides L-type Ca 2+ channels, Na + /Ca + exchanger activity is also found predominantly in the T-tubules of rat ventricular myocytes [23]. 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ current in the T-tubules, since Kawai et al. [24], found Ltype 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 Ca 2+ and Na + ion balance.
Since our study is focused on voltage clamp testing of Ca 2+ 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  Ca 2+ 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 Ca 2+ into the LSR lumen against a concentration gradient. In contrast, the jSR membrane contains an outwardly directed ryanodine (Ry)-sensitive channel for Ca 2+ release from the jSR to the dyadic space. The jSR fluid compartment contains the Ca 2+ 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 Ca 2+ 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 Ca 2+ release channel [25]. 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 Ca 2+ Figure 2 Electrical equivalent circuit for the plasma membrane of a rat ventricular cell. Figure 2: C m, TT : membrane capacitance of the junctional SL membrane coupled with the dyadic space; C m,SL : membrane capacitance of the uncoupled free SL membrane; Currents through the uncoupled free SL membrane are (a) I Ca,L,SL : L-type calcium current, (b) I Na,SL : sodium current through the DHPR channel, (c) I Cs, SL : cesium current through the DHPR channel, (d) I NaCa,SL : sodium-calcium exchanger current, (e) I NaCs,SL : sodium-cesium exchanger current, (f) I PMCA : calcium pump current, (g) I Na,b : background sodium current; Currents through the junctional SL membrane coupled with a dyadic space are (h) I Ca,L,TT : L-type calcium current; (i) I Na,TT : sodium current through the DHPR channel, (j) I Cs,TT : cesium current through the DHPR channel, (k) I NaCa,TT : sodium-calcium exchanger current, (l) I NaCs,TT : sodium-cesium exchanger current; V o : potential in external medium; V i : intracellular potential; The coupling resistance between the surface sarcolemma and the transverse tubules being very small is neglected in our model hence V m : common transmembrane potential across both uncoupled and coupled membranes.
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 Ca 2+ dependent rate functions within the ryanodine receptor model, which affects the open probability P o of the SR Ca 2+ release channel.
With regard to coupling, the DHP and Ry-sensitive Ca 2+ 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 Ca 2+ 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 Ca 2+ channels using Markovian state models ( Figure 4) which include features such as Ca 2+ mediated channel inactivation, a graded CICR process with a "calcium gain" of approximately 6-7, and two-dimensional Ca 2+ diffusion within the dyadic space. Crank [26] 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 Ca 2+ diffusion by solving the 2-D Laplacian   Figure 6). The characteristics of elemental Ca 2+ 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. [27] and Blatter et al. [28]) and duration (full duration at half maximum (FDHM)) which is of the order of 50 ms (reported by Zima et al. [14]). 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 Ca 2+ 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 Ca 2+ from the SR. The measurements of Diaz et al. [29] 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 Ca 2+ release on the particular value of SR Ca 2+ content (i.e., there is a relationship between SR Ca 2+ content  Table 3: Parameters used to model the transmembrane currents I Ca,L , I NaCa , I PMCA , I NaCs and the background sodium current I Na,b . Adopted (*) or estimated ( ‡) from the cited sources. and peak [Ca 2+ ] myo ; Figure 4 Diaz et al. [29]). This plot gives us a glimpse of the input-output relationship of the dyadic coupling unit and indicates that SR Ca 2+ content is an important controlling variable for the CICR process implemented by the DCU model.

L-Type Ca 2+ Current
A multiple state characterization of I Ca,L in rat ventricular myocytes has been reported previously by our group [30]. The gating scheme used in this I Ca,L model has an additional high voltage state C6 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 O2 dhpr . The degree of opening is enhanced in the presence of activated calmodulin-dependent kinase (CaMKII act ) [31][32][33] which is known as Ca 2 + -dependent facilitation (CDF). CaMKII has been shown to tether to the I Ca,L channel [34] functioning as a Ca 2+ signalling sensor for facilitation. The open probability of the I Ca,L channel is also increased in the presence of activated calcineurin (CaN act ) [35].
These two effects are modeled as shown in Figure 4A, where k dhpr 12 is a function of Figure 5 Schematic of the CICR subsystem. Figure 5: Schematic of the CICR subsystem. The dyadic coupling unit comprises of a DHP-sensitive Ca 2+ channel opposing a Ry-sensitive Ca 2+ channel. The transmembrane proteins triadin and junction along with calsequestrin mediate the interaction between the luminal Ca 2+ and the RyR thus regulating the release of Ca 2+ flux from the jSR into the dyadic space.  Figure 4A, state S2 dhpr , which denotes Ca 4 CaM bound to the IQ-motif of the I Ca,L channel, modulates calmodulin dependent Ca 2+ -induced inactivation. This is in agreement with the findings of Nikolai M. Soldatov [37]. 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 [35]. Most of the beat-to-beat modulation is produced by CaM and CaMKII, whereas CaN is constitutively active in the dyad [38]. Tables 5, 6    Our previous study [30] was focused on the characterization of the I Ca,L channel under conditions of low [Ca 2+ ] in the dyadic space and myoplasm. In fact, Ca 2+ -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 Ca 2+ release on I Ca,L , we modified the original I Ca,L model to better characterize the process of Ca 2+ -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 ( k dhpr 25 ), and the Ca 4 CaM dependent inactivation rate function ( k dhpr 24 ). The specific formulas for the Markovian state equations and the modified k dhpr 25 and k dhpr 24 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 Ca 4 CaM-dependent inactivation process. The rate constants k dhpr 36 and k dhpr 63 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 Ca 2+ 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. [39], and Post and Langer [40] considered the effect of Ca 2+ 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) Ca 2+ binding sites on SL wall boundary of cylindrical dyadic space. The presence of these membrane bound sarcolemmal Ca 2+ binding sites in our model has significant physiological implications, in that it prevents local dyadic Ca 2+ concentration near the "mouth" of DHPsensitive 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 Ca 2+ binding sites does not significantly slow the build-up of Ca 2+ within the dyadic space, nor is the Ca 2+ -induced Ca 2+ release mechanism affected [41].

Ca 2+ Release Channel
The gating characteristics of the Ry-sensitive release channel are not only modulated by the dyadic Ca 2+ concentration at its mouth but also the jSR Ca 2+ 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 [42] with 4 Ca 2+ binding for activation and 3 Ca 2+ binding for inactivation to explain the "adaptation" of the RyR observed by Gyorke and Fill [43]; (2) the 6-state scheme suggested by Zahradnikova and Zahradnik [44], which allowed opening the RyR channel upon binding of a single calcium ion; and (3) the Markov model proposed by Keizer and Smith [45], which can be dynamically switched among the six, five and fourstate representations during the simulation as Ca 2+ levels vary. Stern et al. [6] demonstrated that none of these schemes yielded stable local control of SR release, even with extensive adjustment of parameters. Stern et al. [6] 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 [Ca 2+ ] 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 [46] (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 Ca 2+ ions, whereas inactivation depends on the binding of a single Ca 2+ ion. Four additional features have been built into our model: (a) the crucial role of the luminal SR Ca 2+ 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][48][49][50] is modeled via the rate functions k ryr 12 , k ryr 41 , k ryr 43 and k ryr 32 ; (c) a stronger Ca 2+ 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 Ca 2+ concentration consequential to their own Ca 2+ release [51]; and (d) inactivation of the RyR at different depolarization levels is made dependent on local Ca 2+ concentration (i.e., [Ca 2+ ] ryr at the "mouth" of the RyR channel on the dyadic side) as per Wier et al. [52] and Zucchi et al. [53]. Our initial studies indicated that the repriming rate ( k ryr 32 ) in the RyR gating scheme of Stern et al. [46] was quite large, which can lead to a saturated open probability of the Ca 2+ release channel. This occurs during the later phase of channel inactivation, where the large value of k ryr 32 tends to reactivate the channel, resulting in saturated Ca 2+ release. Therefore, we utilized a value of k ryr 32 that is 10% of that used by Stern et al. [46], and further assumed that the unitary permeation flux of the jSR release channel is proportional to jSR luminal Ca 2+ concentration ( [Ca 2+ ] jSR ). The specific equations for the Ca 2+ release model are given in Appendix A1 (Equations 73-92).

Luminal RyR sensor
The RyR Ca 2+ release is modulated by a multi-molecular Ca 2+ signalling complex which is localized to the junctional SR [54][55][56]. This complex consists of the ryanodine    receptor (RyR) which functions as a Ca 2+ conducting pore [57,58], calsequestrin (CS) which acts as the Ca 2+ binding protein [59,60], and the junctional SR transmembrane proteins triadin [61] and junctin [62]. It was known previously that the proteins triadin and junctin anchor calsequestrin to the ryanodine sensitive receptor [61]. More recently, it has been observed that the protein-protein interactions between triadin, calsequestrin and RyR modulate sarcoplasmic reticulum calcium-release in cardiac myocytes [25]. Triadin and junctin are structurally homologous proteins [63], 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. [64] takes a heuristic approach towards free SR Ca 2+ 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. [25] which uncovers complex [Ca 2+ ] jSR ) dependent, CS mediated mechanistic interaction of the protein triadin with the RyR channel. The triadin protein facilitates SR Ca 2+ 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 A1 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 I2 ls . Triadin is also known to exist in a form bound to CS alone [25] denoted by the state I3 ls and in its free unbound form represented as I4 ls in the model. CS which is known to exist in a Ca 2+ bound form modeled by B6 ls also modulates SR Ca 2+ release by influencing the open probability of the RyR channel [65][66][67], via interactions with Triadin [54,66,25]. The degree of unbound CS is denoted by B5 ls in the model.
During the diastolic period the relatively large concentration of available free Ca 2+ in the SR results in most of the CS being bound to Ca 2+ , 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 Ca 2+ induced SR release. Our model incorporates this property by facilitating movement of states towards A1 ls and B6 ls in the presence of large SR Ca 2+ concentration. Following SR Ca 2+ release, a reduced luminal Ca 2+ concentration in the jSR causes an increase in the amount of free calsequestrin (B5 ls in the model) available to bind with triadin. This results in a decrease in the extent of interaction between triadin and RyR (A1 ls in the model), thus inhibiting the RyR channel and leading to robust termination of SR Ca 2+ release in cardiac myocytes [66]. This release termination mechanism is incorporated in our combined RyR-Luminal sensor model ( Figure 4B and Appendix A1, Equations [73][74][75][76][77][78][79][80][81][82][83][84][85][86][87][88][89][90][91][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 A1 ls ) provides adequate peak RyR release which translates into the upstroke velocity of the cytosolic Ca 2+ transient; (ii) the luminal sensor mediated RyR channel inactivation (via A1 ls ) causes timely release termination resulting in a cytosolic Ca 2+ 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 [Ca 2+ ] (both on the dyadic side ( [Ca 2+ ] ryr ) as well as the luminal side ( [Ca 2+ ] jSR )) experiences a strong interaction of Ca 2+ flow through itself and its gating mechanism described as [Ca 2+ ] 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 (A1 ls , I2 ls , I3 ls , I4 ls ) and CS (I2 ls , I3 ls , B5 ls , B6 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.

Ca 2+ Buffering in Myoplasm and SR
Ca 2+ buffers play an important role in sequestering a fraction of the total Ca 2+ 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 Ca 2+ binding to calmodulin were based on a model from [22], whereas those for Ca 2+ binding to troponin were taken from Potter and Zott [72]. Rate constants used to describe Ca 2+ binding to calsequestrin were based on the study of Cannell and Allen [73], whereas those for Ca 2+ binding to troponin-Mg complex were adopted from Lindblad et al. [74].
It is well recognized that fluorescent indicator dyes introduced into the cytosol also act as Ca 2+ buffers [22], even at submillimolar concentrations [51]. 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. [75].

Ca 2+ -Uptake
Cytosolic Ca 2+ is pumped into the LSR ( Figure 1A) by a Ca 2+ -ATPase, which lowers [Ca 2+ ] myo and helps to induce relaxation in cardiac muscle. The transport reaction involves two Ca 2+ ions and one ATP molecule [76], 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. [77] and takes into account both the forward flux of Ca 2+ from the cytosol to the LSR lumen and the backward flux from the LSR to the cytosol along with the Ca 2+ 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 Ca 2+ transport by increasing the pumping rate [78] and indirectly via phosphorylation of PLB [79] relieving the inhibition caused by PLB on the SERCA pump in turn increasing the sensitivity of the pump for Ca 2+ uptake. These two effects are modeled by allowing the rate constants for Ca 2+ + 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 Ca 2+ between the cytosol and the SERCA protein. Similarly, the current I serca,sr dictates the transport of Ca 2+ between the SERCA protein and the LSR. The difference in these Ca 2+ currents accounts for the Ca 2+ buffered by the SERCA protein. The jSR is subsequently refilled by Ca 2+ diffusion from the LSR. The differential equations describing the Ca 2+ balance and particularly the transfer between two separate SR compartments (jSR and LSR) are provided in Appendix A1 (Eq. 72).

Ca 2+ -Extrusion via Sarcolemmal Ca 2+ Pump
Although the sarcolemmal Ca 2+ -pump has a high affinity for [Ca 2+ ] myo , its transport rate is far too slow for it to be an important factor in Ca 2+ fluxes during the cardiac cycle. It might, however be more important in long-term extrusion of Ca 2+ by the cell. Our model of the plasma membrane Ca 2+ ATPase pump current is adopted from Sun et al. [30]. 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 [82].

Ca 2+ -Extrusion via Na + /Ca 2+ Exchanger
In mammalian cardiac cells, it is generally accepted that the Na + /Ca 2+ exchanger has a stoichiometry of 3Na + :1Ca 2+ [83]. I NaCa (again, the combination of I NaCa,TT and I NaCa, SL ) is important in removing Ca 2+ during twitch relaxation, in competition with I up . A simple thermodynamic Na + /Ca 2+ exchanger current model [84] may be sufficient to predict the direction of Ca 2+ transport by Na + /Ca 2+ exchange and the driving force, however the amplitude is subject to kinetic limitations (depending on substrate concentrations). A more comprehensive Na + /Ca 2+ exchanger current equation [85] is adopted in our model, which includes factors for allosteric Ca 2+ activation and the transport for Na + and Ca 2+ inside and out. The maximal flux through the exchanger V max is estimated to ensure that the Ca 2+ ion transport (which is voltage dependent) via the Na + /Ca 2+ exchanger matches the influx of Ca 2+ via the I Ca,L channel [86,87], maintaining whole cell Ca 2+ homeostasis.

Na + /Cs + Pump
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 [88]. However, in our experimental protocol, external solution in the bath was normal Tyrode (1 mM Ca 2+ ) 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 Na 2 ATP , 3.5 mM MgCl 2 and 5 mM HEPES. Activation of the electrogenic sodium pump in mammalian non-myelinated nerve fibres [89], skeletal muscle [90] and rat brain cells [91] 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 [92]. Monovalent cations (including Cs + ) were also added to K + free bathing solution to reactivate the sodium pump in guinea-pig ventricular myocytes [93]. The effectiveness of Cs + as an external cation in activating the electrogenic sodium pump is known to be lesser than potassium [92].
We have represented I NaCs in our model by the expression for Na + /K + pump formulated by Linblad et al. [74] 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. [94].

Results
The DCU is a fundamental element in the mechanism of CICR. The sequence of events resulting in CICR is triggered by the Ca 2+ entering through the I Ca,L channel. The characteristics of this channel are hence examined in detail to understand its voltage and Ca 2+ 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 Ca 2+ mediated interaction between these channels is investigated in detail. This is followed by an effort to understand the role of [Ca 2+ ] jSR in CICR. In particular, the relationship of the peak cytosolic calcium transient and the SR Ca 2+ 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 Ca 2+ current (I Ca,L )
The L-type DHP-sensitive Ca 2+ 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 Ca 2+ concentration in the vicinity of the channel [95]. The relative contribution of the voltage-dependent and Ca 2+ -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 Ca 2+dependent inactivation via I Ca,L self-inhibition in the absence of RyR release; and (iii) Ca 2+ substituted with Ba 2+ to facilitate only the voltage-dependent inactivation pathway and block all Ca 2+ dependent inactivation pathways (k24 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 [18]). The degree of I Ca,L channel opening is known to be enhanced in the presence of CaMKII act [33] which is known as Ca 2+ -dependent facilitation (CDF). According to Bers et al. [33] the positive regulation of the channel by CaMKII act requires Ca 2+ influx (it is not seen when Ba 2+ is the charge carrier as in Figure 7D and is more strongly apparent when local Ca 2+ influx is amplified by SR Ca 2+ release as in Figure 7A). Increase in CaMKII act seems to play an important function in enhancing I Ca,L peak amplitude [96], suggesting a critical role for Ca 2+ -dependent facilitation. Our model generated results also indicate that a stronger CDF (in Figure 7A) in the presence of RyR Ca 2+ 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 (Ca 2+ substituted by Ba 2+ )). 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, Ca 2+ -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, Ca 2+ -dependent inactivation (CDI) caused by Ca 2+ release from Ry-sensitive channels is much more significant than the self-inhibition produced by Ca 2+ influx via I Ca,L itself. The graded behavior of the cytosolic Ca 2+ 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 Ca 2+ inactivation effects have been blocked by Ba 2+ substitution for Ca 2+ . 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 O2 dhpr -C4 dhpr pathway, owing to the large Ca 2+ 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 O2 dhpr -C4 dhpr pathway, because of the low cytoplasmic Ca 2+ 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 Ca 2+ 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 [18]) 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 Ca 2+ transient is shown in Figure 8C. Figure  8B shows the well known [52] bell-shaped dependence of the peak [Ca 2+ ] 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 Ca 2+ 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; Ca 2+ entering via the I Ca,L channel causes Ca 2+ induced inactivation, although the magnitude of inhibition is far less than in the control case. Case 3: Ba 2+ substitution for Ca 2+ to completely suppress the Ca 2+ inactivation mechanism; the only inactivation pathway present is the slow voltage dependent pathway [97], 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 Ca 2+ release (a result of enhanced CDF in the presence of elevated Ca 2+ levels). This highlights the important role of Ca 2+ in regulating I Ca,L channel opening and thus controlling the amount of extracellular Ca 2+ entering the cell.
One of the most important characteristics of the CICR mechanism is the modulation of the RyR Ca 2+ 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 Ca 2+ release. As shown in Figure 10A, our model reproduces graded release. It is important to note that the onset of the Ca 2+ 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 Ca 2+ release channel is controlled by the amount of trigger Ca 2+ present at the 'mouth' of this channel, which in-turn is graded by the amount of trigger Ca 2+ 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 Ca 2+ release channel ( Figure  4) to be [Ca 2+ ] ryr dependent. Small changes in i Ca,L,TT cause modulation of pre-release [Ca 2+ ] ryr , which dictates the propensity of the RyR channel for a Ca 2+ release. The pre-release values of the rate constants in the 4-state Markovian model of the RyR channel (nonlinear dependence on [Ca 2+ ] ryr as shown in Appendix A1, Equations 79-84) along with the diastolic [Ca 2+ ] jSR (which is kept constant in Figure 10) set the peak RyR open probability achieved by the RyR channel as well as the peak [Ca 2+ ] myo . The use-dependent adaptation of the RyR channel [13] is reflected in its non-linear response to trigger Ca 2+ . As shown in Figures 10B-iii, iv, the delay in the onset of cytosolic Ca 2+ transient closely follows the modulation of the onset of RyR release. Though occurring in separate compartments, the peak of cytosolic Ca 2+ transient also tracks the maximum value attained by the open probability of the Ry-sensitive Ca 2+ release channel. Here, the clamp voltage was held constant at 10 mv to avoid the interference of voltage dependent Ca 2+ transport via the Na + /Ca 2+ exchanger which is co-located [98] 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 Ca 2+ 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. [13] bringing trigger Ca 2+ into a unitary dyad. This L-type Ca 2+ channel (LCC) triggers release from a cluster of RyR channels causing a Ca 2+ spark ( Figure 6), which results in the rapid increase in [Ca 2+ ] 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 Ca 2+ spark increases with decrease in the magnitude of i Ca,L,TT [13]. The SR Ca 2+ release flux underlying a typical Ca 2+ spark corresponds to approximately 2 pA [27,28,13]. From RyR single channel conductance measurements in lipid bilayer studies [99], a Ca 2+ spark translates into a release from around 4 RyR channels (also reported by Blatter et. al. [28]). This single Ca 2+ 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 Ca 2+ release flux into the cytosol causing the Ca 2+ 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 [25] that the luminal sensor plays a fundamental role in an active extinguishing mechanism [51] that effects a robust [Ca 2+ ] jSR -dependent closure of the RyR channel. This mechanism for inactivation of the Rysensitive Ca 2+ release channel is accounted for in our model. Cytosolic Ca 2+ 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][101][102]. The effect of depolarization duration on the time course of the cytosolic Ca 2+ 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 Ca 2+ transient. This effect is a combined result of (i) the modulation of release due to pulse duration dependent change in the amount of trigger Ca 2+ entering the dyadic coupling unit (DCU) and (ii) the pulse duration dependent change in the relative role of the Na + /Ca 2+ exchanger co-located in the DCU.

High Gain of Ca 2+ 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 Ca 2+ entering the DCU via the I Ca,L,TT channel causes a large release of Ca 2+ from the sarcoplasmic stores via the Rysensitive Ca 2+ release channel on the jSR lumen that interacts with the DCU. This high gain is essential in producing physiological cytosolic Ca 2+ 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: 1. average integrated RyR flux to average integrated DHPR flux [7]; and 2. the peak Ca 2+ transient in the presence of CICR to the peak calcium transient in its absence, con-tributed by the trigger calcium alone [3].
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 Ca 2+ flux through DHPR and RyRs, respectively. Figure 12B shows cytosolic calcium transients under conditions of ryanodine blockade and with RyR Ca 2+ 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 [3]. This result is also consistent with gain calculations from the measured data of Fan and Palade [17] on rat ventricular cells. They estimated a gain of approximately 7 using comparisons of the rates of rise of Ca 2+ 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. [104]). It is important to note that, with increasing clamp voltage, the decreased ability of the Na + /Ca 2+ exchanger (which is co-located [98] in the dyad) to extrude Ca 2+ 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 [105]. It is now known [25] that the luminal sensor plays a fundamental role in an active extinguishing mechanism [51] effecting a [Ca 2+ ] 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, [Ca 2+ ] ryr drops to a level much below what is caused by a sparklet (trigger Ca 2+ ) 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 Ca 2+ 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 [Ca 2+ ] 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 (T 2 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 (T 1 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 (T 2 ) 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 [Ca 2+ ] 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 Ca 2+ release based on the degree of SR Ca 2+ content recovery. This is an important feature that helps restore the jSR Ca 2+ 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 [106] on rat ventricular myocytes. When the RyR channel is in the absolute refractory period, [Ca 2+ ] RyR drops to a level much below what is caused by a sparklet (trigger Ca 2+ ) 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 [Ca 2+ ] myo and [Ca 2+ ] 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 Ca 2+ 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 Ca 2+ concentration consequential to its own  Ca 2+ release [51]. The lack of a luminal sensor forces the RyR channel to solely rely on Ca 2+ dependent inactivation mechanism. The resulting inadequate RyR inactivation depletes diastolic [Ca 2+ ] jSR level ( Figure 14B-ii) causing the cytosolic Ca 2+ 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 [Ca 2+ ] 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 Ca 2+ content, is key to long term Ca 2+ stability in the cell.

CICR modulation by the jSR Ca 2+ content
From the refractory characteristics of the RyR channel, it is evident that the RyR release is strongly dependent on the jSR Ca 2+ content. In fact, the SR Ca   Figure 7), hence causing an indirect self-inhibition of RyR release. Figure 15 shows a phase plot of RyR open probability (O2 ryr ) versus the Ca 2+ concentration at the mouth of the RyR channel in the dyadic space [Ca 2+ ] 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 [Ca 2+ ] jSR is 918 μM. Following excitation by the trigger current during phase A, the RyR channel begins to open, allowing Ca 2+ flux into the dyad and elevating [Ca 2+ ] ryr at the mouth of the RyR channel. The rate of opening of the channel is Ca 2+ dependent; hence as [Ca 2+ ] ryr increases, RyR open probability increases first slowly and then ever more rapidly. With this onrush of Ca 2+ , [Ca 2+ ] ryr rapidly equilibrates with the [Ca 2+ ] jSR , reducing and eventually abolishing the concentration gradient across the channel. [Ca 2+ ] ryr soon reaches its maximum value (T1 in Figure 15) and begins to track the decrease in [Ca 2+ ] jSR , there being no further significant concentration gradient existing between [Ca 2+ ] jSR and [Ca 2+ ] ryr at the mouth of the RyR channel. During phase B, the rate of rise in RyR open probability decreases due to decreasing [Ca 2+ ] ryr and the onset of self-inhibition due to the large values of Ca 2+ 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 [Ca 2+ ] dyad levels as a result of (a) Ca 2+ fluxing out of the dyad into the cytosol, (b) lack of a drive from SR release due to drastically diminished Ca 2+ gradient Figure 15 Ryanodine receptor open probability. Figure 15: RyR open probability modulation by concentration at the mouth. The pre-release diastolic [Ca 2+ ] jSR is 918 μM. The stimulation protocol used is a pulse of amplitude (-40 mv to 10 mv) and duration (50 ms). Phase A -As a result of excitation by the trigger current, [Ca 2+ ] ryr initially increases rapidly with a small accompanying rise in RyR open probability followed by steeper increase to peak [Ca 2+ ] ryr with a large increase in RyR open probability.; Phase B -RyR open probability increases to its maximum value at a declining rate reflecting a falling [Ca 2+ ] ryr level and onset of Ca 2+ induced self-inhibition; Phase C -Decreasing Ca 2+ dyad and continued Ca 2+ induced selfinhibition, begins the channel recovery process; Phase D -Decreasing [Ca 2+ ] jSR forces a strong RyR inactivation via the luminal sensor hence robustly closing the channel. at the mouth of the RyR channel and (c) continued Ca 2+ induced self-inhibition. Decreasing levels of free Ca 2+ 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 [Ca 2+ ] jSR occurs at T5, and the channel is ultimately closed at T6 (Figure 15). The luminal sensor mediated inactivation assists in robust [Ca 2+ ] jSR recovery.
Cytosolic peak [Ca 2+ ] dependence on SR Ca 2+ content Figure 16A shows three characteristic regions in the plot of the peak [Ca 2+ ] myo vs [Ca 2+ ] jSR . In Region I, where the jSR load is small, two things occur: (a) RyR release is reduced, as is the basal intracellular Ca 2+ concentration, and (b) Ca 2+ dependent inactivation of the I Ca,L,TT channel is reduced due to the lowered dyadic Ca 2+ concentration, which allows for greater Ca 2+ entry into the dyadic space via the I Ca,L channel. This increase in trigger current results only in a very small increase in Ca 2+ release, because of the inhibition of the Ry-sensitive channel caused by the luminal sensor in response to low SR Ca 2+ content. Thus, these factors cumulatively result in only a small linear increase in the peak [Ca 2+ ] myo with increasing SR Ca 2+ content in Region I. The behavior in Region I corresponds to traces 7-1 in Figure 16B (shows phase plots of O2 ryr vs [Ca 2+ ] ryr for different steady state [Ca 2+ ] jSR concentrations), where the peak RyR open probability shows very small gradual increase with increasing levels of diastolic pre-release [Ca 2+ ] jSR for a constant trigger Ca 2+ input via I Ca,L . The steady state diastolic [Ca 2+ ] jSR levels being low, the maximum RyR open probability values attained are low. Although the peak [Ca 2+ ] ryr does not increase substantially as the diastolic [Ca 2+ ] jSR levels increase from 317μM to 673μM, the area enclosed by the [Ca 2+ ] ryr -O2 ryr loop (which indicates the amount of SR Ca 2+ 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 Ca 2+ content begins to translate into a substantial increase in peak [Ca 2+ ] myo . This is because increasing diastolic [Ca 2+ ] 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 Ca 2+ transient. The inhibiting role of the luminal sensor is relieved with building [Ca 2+ ] jSR levels. The nonlinearity observed between the diastolic [Ca 2+ ] 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 [Ca 2+ ] jSR levels increase beyond 734μM, the area enclosed by the [Ca 2+ ] ryr -O2 ryr loop increases substantially reaching a maximum value at [Ca 2+ ] jSR level of 948μM not only due to a significant increase in both peak [Ca 2+ ] ryr and peak O2 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 [Ca 2+ ] myo due to the large RyR Ca 2+ release associated with the high SR Ca 2+ 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 [Ca 2+ ] jSR values larger than 850μM (traces 6-1 in Figure 16C). The continuing rapid increase in peak cytosolic Ca 2+ concentration in region III is a combined result of: (a) increased release owing to large SR Ca 2+ 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 Ca 2+ concentration levels in the cytosol after RyR release [77]. As the diastolic [Ca 2+ ] jSR levels increase beyond 948μM, the area enclosed by the [Ca 2+ ] ryr -O2 ryr loop begins to decrease despite increasing peak [Ca 2+ ] ryr due to saturation of the peak RyR open probability and increasing rate of RyR channel inactivation due to large values of local Ca 2+ concentration ( [Ca 2+ ] ryr ) assisting in faster recovery [51]. This model generated relationship between the peak [Ca 2+ ] myo and the pre-release [Ca 2+ ] jSR agrees with Figure 4 in Diaz et. al. [29].
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 [Ca 2+ ] ryr the RyR open probability takes a lower value despite a higher steady state diastolic [Ca 2+ ] jSR level. This decreasing slope with increasing diastolic, pre-release [Ca 2+ ] jSR reflects a faster rate of rise in [Ca 2+ ] ryr for a larger pre-release SR content with a very small accompanying increase in the RyR open probability.

Caffeine
Adachi-Akahane et al. [71] investigated I Ca,L -induced Ca 2+ release in rat ventricular myocytes in the presence and absence of caffeine, which altered [Ca 2+ ] myo . Using their protocols and data as a guide, we modeled the effect of caffeine on jSR Ca 2+ 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 Ca 2+ release. The process is assumed to be dosedependent, 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 [Ca 2+ ] 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 [Ca 2+ ] 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 Ca 2+ 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 [Ca 2+ ] 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 Figure 17 Effects of the application of caffeine. Figure 17: Model-simulated effects of the application of caffeine on I Ca,L,TT and [Ca 2+ ] i in a cell dialyzed with 2 mM Fura-2. Changes in [Ca 2+ ] myo resulting from: (A) application of a simple depolarizing pulse from -40 mV to 0 mV, and (B) the same depolarizing pulse following pre-exposure to caffeine (5 mM). Panel (C) compares the I Ca,L waveforms associated with these protocols. This model-generated result shows similarity to data in Figure 3, Adachi-Akahane et al. [71].
smaller [Ca 2+ ] myo transient via CICR (30 nM, Figure 17B, second peak). This I Ca,Levoked Ca 2+ 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 [Ca 2+ ] jSR ([Ca 2+ ] 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 [Ca 2+ ] dhp during release does not increase as much as in the control case, since jSR Ca 2+ release is reduced. Consequently, Ca 2+ 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 Ca 2+ concentration at the mouth of the DHP channel ( [Ca 2+ ] 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. [41], that Ca 2+ induced inactivation is a result of calmodulin mediated sensing of the local Ca 2+ concentration.

Thapsigargin
Analogous to the experimental protocol of Adachi-Akahane et al. [71] 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 Ca 2+ uptake is blocked, SR Ca 2+ content declines, as does the jSR Ca 2+ release with each voltage pulse ( Figure 18A). Considerable time is consumed after thapsigargin application, before SR Ca 2+ content is finally exhausted (8-10 beats at 0.05 Hz; Figure 18A). As the magnitude of the [Ca 2+ ] myo transient decreases, the [Ca 2+ ] dhp dependent inactivation of I Ca,L slows ( Figure 18B), resulting in an enhanced peak.
For a given clamp amplitude, a decrease in jSR Ca 2+ content decreases Ca 2+ release, producing decreased Ca 2+ concentration at the dyadic side of the DHP-sensitive channel. Consequently, a smaller amount of calcium-calmodulin complex (Ca 4 CaM) is activated, leading to a lower degree of Ca 2+ 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 Ca 2+ entering the cell. This inverse relationship between jSR Ca 2+ content and the amount of trigger Ca 2+ 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 Ca 2+ -buffers, since the small dimensions of the dyadic space allow jSR released Ca 2+ to reach its destination (i.e., the DHP receptor) before being captured by the buffers [109].

Effect of modulation of [Ca 2+ ] o
Increase in [Ca 2+ ] o causes an increase in entry of extracellular Ca 2+ into the cell via the I Ca,L channel. This results in an increase in SR Ca 2+ content, as shown in Figure  19B. Aided by an increase in trigger Ca 2+ as well as a larger SR Ca 2+ content, the RyR release is enhanced. This increase in SR release assisted by impaired Na + /Ca 2+ exchange due to elevated [Ca 2+ ] o manifests in an increase in the peak Ca 2+ concentration levels in the cytosol, as shown in Figure 19A. Similarly, a decrease in extracellular Ca 2+ levels causes a decrease in trigger Ca 2+ available to cause CICR, which results in a decrease in peak [Ca 2+ ] levels in the cytosol, as shown in Figure 19C. Reduced Ca 2+ entry into the cell via I Ca,L , translating into decreased availability of post-release Ca 2+ in the cytosol, causes reduction in SR filling via the bi-directional SERCA pump. The corresponding drop in the SR Ca 2+ content is evident in Figure 19D. Although, our model-generated results show similarity to data in Figure 2, Diaz et al. [110], our result  This model-generated result shows similarity to data in Figure 11, Adachi-Akahane et al. [71] in Figure 19D does not confirm an increase in SR Ca 2+ content with a decrease in [Ca 2+ ] o as reported in Figure 2B, Diaz et al. [110]. 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 Ca 2+ content with decrease in [Ca 2+ ] o in cells not showing spontaneous release. In this study, we have attempted to isolate the effects of changes in [Ca 2+ ] o on intracellular Ca 2+ levels by ensuring that there are no voltage dependent effects on Ca 2+ transport via the I Ca,L channel, Na + /Ca 2+ exchanger or the plasma membrane Ca 2+ ATPase pump. The simulation protocol used in our study (reproduced from Diaz et al. [110]), 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 Ca 2+ concentration, in response to a voltage clamp pulse train, increasing clamp potential (>10 mv) results in a decrease in steady state SR Ca 2+ content as a result of the dominating effect of reduced Ca 2+ entry via the I Ca,L channel despite a decrease in net Ca 2+ extrusion via the Na + /Ca 2+ exchanger per cycle.

Calcium balance under conditions of repetitive stimulation
Our model exhibits long term Ca 2+ stability. Analogous to the experimental protocols of Negretti el al. [111], 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 Figure 19 Effects of modulating [Ca 2+ ] o . Figure 19: The effects of changing external Ca 2+ concentration on cytoplasmic Ca 2+ and SR Ca 2+ content. The stimulation protocol used is a pulse train of amplitude (-40 mv to 0 mv), duration (100 ms) and frequency of 0.5 Hz. (A) Following sudden rise in [Ca 2+ ] o the instantaneous increase in [Ca 2+ ] myo occurs due to more Ca entry, followed by a gradual staircase rise due to the the change in SR content; (B) Changes in [Ca 2+ ] myo forced by a sudden decrease in [Ca 2+ ] o ; (C) Changes in [Ca 2+ ] jSR due to increased [Ca 2+ ] o ; (D) Changes in [Ca 2+ ] jSR due to reduced [Ca 2+ ] o . This modelgenerated result shows similarity to data in Figure 2, Diaz et al. [110]. steady-state operation, Ca 2+ entry into the cytosol via I Ca,L must exactly balance Ca 2+ efflux. Changing the rate or pattern of stimulation can have significant effects on the cell's cytosolic Ca 2+ balance and subsequent contractile response [16]. To demonstrate Ca 2+ 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 [Ca 2+ ] myo were correlated with sarcolemmal and SR Ca 2+ fluxes. Analogous to the protocol of Negretti et al. [111], 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 (shortpulse) at a frequency of 0.33 Hz, for a total test duration of 2 minutes (40 pulses). In analyzing the Ca 2+ balance, the amount of electric charge Q transferred per pulse was Numerical integration was carried out with respect to steady-state levels, hence only activating currents and leakage currents were considered (background Ca 2+ current is neglected). The results of the Ca 2+ balance are presented in a manner where relative Ca 2+ 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 Ca 2+ (the exchanger

Long-pulse protocol
During repetitive long-pulse (800 ms) stimulation, transient Ca 2+ fluxes cross the sarcolemmal and SR membranes in both directions (inward (into the cytoplasm) and outward). However, [Ca 2+ ] 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 Ca 2+ influx and efflux to and from the cytoplasm, under control conditions. In contrast, this Ca 2+ balance is not present in the shortpulse protocol, and [Ca 2+ ] myo is elevated at pulse termination. Our long-pulse simulations show that: (1) the magnitude of individual Ca 2+ 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 Ca 2+ buffer in the jSR compartment is reduced by 15% as a result of Ca 2+ release during each cycle. Although Ca 2+ 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 Ca 2+ nor consecutive-pulse Ca 2 + depletion (or augmentation) in the jSR compartment under the long-pulse protocol.
In Figure 20, we refer to Ca 2+ 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 Ca 2+ 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 Ca 2+ charge entering the cell via I Ca,L is 603.28 pC (pico-coulomb), and the sum of averaged Ca 2+ 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 Ca 2+ 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 Ca 2+ via the trigger current I Ca,L for a prolonged duration of 800 ms compared to 50 ms previously. The integral of SR Ca 2+ uptake is 751.8 pC. These two internal component currents yield flat integral values during long-pulse stimulations and sum to zero.

Short-pulse protocol
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 [Ca 2+ ] myo observed with stimulation by long pulses, short pulse stimulation produces a negative staircase ( Figure 20B) in agreement with the data (Negretti et al. [111]), 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 Ca 2+ concentration [Ca 2+ ] jSR . Both of these indicators show a net loss of Ca 2+ from the cell. Figure 20C and 20D show the transient charge movements associated with the external and internal Ca 2+ currents respectively during application of the short-pulse train, which are indicative of an elevated Ca 2+ load and resultant Ca 2+ imbalance when the pulse train is first applied. Balance is achieved, but at new lower values of [Ca 2+ ] myo and [Ca 2+ ] 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 Ca 2+ load on the cell. The decreasing amplitude of the Ca 2+ transient causes a decreased LSR uptake current (I cyt,serca ), which in turn causes a decrease in [Ca 2+ ] 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 Ca 2+ 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 Ca 2+ levels of [Ca 2+ ] myo over this prolonged time frame are indicative of long tern Ca 2+ stability in the model. The corresponding [Ca 2+ ] 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 [Ca 2+ ] 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 Ca 2+ influx called a "tail current", which has been shown to trigger RyR release and hence cause contraction during repolarization [100,2]. The secondary Ca 2+ 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 Ca 2+ release [112]. Our model reproduces this secondary Ca 2+ 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. ([8]; 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 Ca 2+ channel, with RyR Ca 2+ release indicated by the corresponding decrease in [Ca 2+ ] jSR (panel C) and increase in [Ca 2+ ] myo (panel D). We note from Figure 21C that for a shorter 50 ms VC pulse, there is no RyR Ca 2+ release although an I Ca,L tail current is produced. This tail current produced with the 50 ms pulse slightly augments the [Ca 2+ ] myo transient just after its peak (trace 1, panel E), but there is no corresponding decrease in the [Ca 2+ ] jSR transient (trace 1, panel D) indicating no RyR Ca 2+ 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 [Ca 2+ ] myo transient as seen in data corresponding to trace 4 of Figure 7C is a result of the large Ca 2+ influx via the I Ca,L "tail current". During this phase, the Na + /Ca 2+ exchanger is biased to extrude Ca 2+ out of the dyad, and hence cannot be responsible for this secondary rise in the cytosolic Ca 2+ transient. As pulse duration is increased beyond 100 ms the "tail current" causes a

P4
Fraction of active Thr 287 -autophosphorylated but Ca4CaM trapped CaMKII 9.121920 × 10 -8 Table 10: Initial conditions used for the state vector in order to solve the linearized system of ordinary differential equations.  Figure 21) indicating progressive recovery from the refractory period. Minimum value of [Ca 2+ ] jSR in traces 3-7 of Figure 21D show increase in SR Ca 2+ release with increasing pulse duration as a result of adequate RyR channel recovery and increasing I Ca,L "tail current". Cytosolic Ca 2+ transients shown in Figure 21E indicate "tail current" induced secondary RyR release for pulse duration ≥ 100 ms. It has also been previously reported [113] that recruitment of an additional population of previously 'silent' Ca 2+ 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 O2 dhpr &O3 dhpr ) to the tail current during repolarization but also allows for the voltage-and time-dependent recruitment of 'silent' Ca 2+ channels modeled using a high voltage state C6 dhpr . This study proves to be an adjunct to the study presented in Figure 13 in establishing the refractory nature of the RyR channel.

Cytosolic Buffering
In our model, cytosolic buffering is attributed to several factors including: (a) calmodulin (CaM); (b) the Ca 2+ -specific (Tc) troponin binding site; (c) the Ca 2+ -Mg 2+ competitive troponin binding site; and (d) the fluorescent indicator dye Fluo3 used to detect changes in [Ca 2+ ] myo . The effects of the component buffers in helping to maintain a low [Ca 2+ ] myo are shown in Figure 22 in terms of occupancy functions, such as O C Figure 22 Occupancy of Ca 2+ buffers. Figure 22:  Figure 22 shows that when the Ry-sensitive Ca 2+ release channel is triggered, the jSR releases Ca 2+ and the occupancy O calse in the jSR, declines from 68% to 37%. Ca 2+ release from the jSR induces fast Ca 2+ 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 Ca 2+ and the troponin-Mg sites result in an increase in occupancy of these sites by Ca 2+ (O tmgc ), and a decrease in occupancy by Mg 2+ (O tmgmg ).

Discussion
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 Ca 2+ 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 Ca 2+ 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 Ca 2+ channels separated by a small dyadic space ( Figure 1B): the sarcolemmal DHPsensitive "trigger" channel and the Ry-sensitive jSR Ca 2+ release channel. The trigger channel is voltage-activated and is driven in our simulations by a voltage clamp pulse which opens the channels, admitting Ca 2+ influx (trigger current) to the dyadic space. After diffusion within the dyadic space, trigger Ca 2+ affects the Ca 2+ dependent open probability of the jSR Ca 2+ release channel over a small range of Ca 2+ concentrations. Clamp voltage magnitude strongly affects the I Ca,L current and hence the amount of Ca 2+ delivered to the dyadic space, but the range over which the Ca 2+ concentration at the mouth of the Ca 2+ 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 Ca 2+ 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 Ca 2+ 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 Ca 2+ . Hence the luminal sensor, which is a key element responsible for robust post-release RyR inactivation and refractoriness of the Ry-sensitive Ca 2+ release channel, is critical in providing realistic fits to cytosolic Ca 2+ transients and an adequate refilling time for the SR Ca 2+ stores. Ca 2+ release is thus graded with Ca 2+ concentration at the mouth of Ca 2+ 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 Ca 2+ 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 Ca 2+ release in CICR. It contains Ry-sensitive Ca 2+ channels that require a Ca 2+ concentration gradient directed across the channel and into the dyadic space for their operation. Thus, jSR Ca 2+ concentration must be maintained within an acceptable range so that calcium is always available for ready release. The relationship between jSR Ca 2+ content and peak [Ca 2+ ] 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 Ca 2+ from the myoplasm. Pumping rate is controlled by the cytosolic Ca 2+ concentration [Ca 2+ ] myo as well as the LSR Ca 2+ concentration [Ca 2+ ] LSR as part of a mechanism for replenishing and maintaining jSR Ca 2+ stores. Our model also incorporates the Ca 2+ induced CaM mediated effects of CaMKII and CaN on targets such as the DHP-sensitive I Ca,L channel, Ry-sensitive Ca 2+ release channel, as well as the SR Ca 2+ 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 Ca 2+ signaling. Although physically separated, the voltage-sensitive Ca 2+ channel ( Figure 9) as well as the Ry-sensitive release channel (Phase C in Figure 15) are inhibited by increased Ca 2+ levels in the dyadic space [53]. In earlier models which lacked the luminal sensor, RyR self-inhibition by dyadic Ca 2+ 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 Ca 2+ concentration. Our model predicts (data not shown) that this Ca 2+ signaling continues even in the presence of high concentrations of Ca 2+ buffering agents in the cytosol (in agreement with the data of Adachi-Akahane et al. [71] and Diaz et al. [110]). Tight coupling between the DHP and Ry-sensitive Ca 2+ channels within the dyadic space thus preserves the mechanism of CICR under these extreme conditions.
Regulation of cytosolic Ca 2+ concentration [Ca 2+ ] myo is evident at the cellular level. It is affected strongly by the voltage and [Ca 2+ ] myo -dependent properties of the sarcolemmal currents I NaCa and I Ca,L , as well as the [Ca 2+ ] myo -dependent properties of the I PMCA and SERCA pumps. We demonstrate a model-generated whole-cell Ca 2+ balance, which shows the importance of the Na + /Ca + exchanger in extruding the Ca 2+ that has entered the cell under normal activity, and also any excess that might occur when cytosolic Ca 2+ levels rise. The variety of experiments emulated in this study demonstrates quantification of Ca 2+ balances for all external and internal Ca 2+ fluxes and shows that the model has long-term stability in regulating cytosolic Ca 2+ , 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 Ca 2+ -dependent inactivation as a function of jSR Ca 2+ content ( Figure 17C). Trigger current is voltage-dependent and therefore parameterized by clamp pulse amplitude (Figure 7A). In addition, its inactivation is Ca 2+ dependent, especially during Ca 2+ release. At a constant depolarizing voltage pulse, the pre-release jSR Ca 2+ content dictates the peak of the cytosolic Ca 2+ transient ( Figure 16A), by controlling peak Ca 2+ concentration in individual dyads. As jSR Ca 2+ content increases, so does the peak dyadic Ca 2+ concentration, and the amount of sarcolemmal Ca 2+ influx declines due to greater Ca 2+ dependent inactivation. This autoregulatory feedback mechanism helps to establish a stable operating point for jSR Ca 2+ content. Transient increases in jSR Ca 2+ content bring about increases in Ca 2+ release, but reflexly, decrease sarcolemmal Ca 2+ influx via I Ca,L . Coupled with other dyadic and extra-dyadic mechanisms (I NaCa ) that decrease [Ca 2+ ] myo and hence Ca 2+ -uptake via the SERCA pump, jSR Ca 2+ content decreases. For decreases in jSR Ca 2+ content, the opposite occurs, so that regardless of the sign of the perturbation in jSR Ca 2+ 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 Ca 2+ transient vs SR Ca 2+ 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 Ca 2+ content (simulated either by a decrease of Ca 2+ uptake into the LSR by thapsigargin or an increase of Ca 2+ leak out of the jSR by caffeine) decreases systolic [Ca 2+ ] 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 Ca 2+ transients are commonly observed [114].

Model Limitations
(a) This model of a rat ventricular myocyte is limited to Ca 2+ 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 Ca 2+ dynamics allows one to comprehend more clearly the important role of Ca 2+ signalling pathways and feedback control systems in maintaining whole cell homeostasis over a prolonged period of time. (b) A single dyadic space in our model has one representative, lumped DHP-sensitive and Ry-sensitive Ca 2+ 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 Ca 2+ 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 Ca 2+ transients. The effectiveness of this model is further demonstrated by its ability to accurately characterize the interaction between the DHP and Ry-sensitive Ca 2+ channels, including pulse duration dependent termination of release (Figure 11), Ca 2 + dependent inactivation of I Ca,L (Figures 7, 9), as well as the wide variety of wholecell voltage clamp protocols (Figures 17, 18, 19 and 20). (c) Although our model provides secondary Ca 2+ 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. ([8]; Figure 8). The restitution time for the RyR channel in rat ventricular myocytes is believed to be at least 25 ms [117] and as fast as 150 ms [118] 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 Ca 2+ spark restitution in rat ventricular myocytes reported by Sobie et al. [106]. 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.

Conclusion
We have developed a mathematical model of Ca 2+ 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 Ca 2+ 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 Ca 2+ release channel. This element is critical in providing realistic fits to cytosolic Ca 2+ transients and an adequate refilling time for the SR Ca 2+ stores. Our voltage-clamp simulations demonstrate graded Ca 2+ transients with sufficient gain, as well as quantification of Ca 2+ balances for all external and internal Ca 2+ 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 Ca 2+ , 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 Ca 2+ 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.

Appendix
Below is the complete set of equations used in the model. . (6) ks dhpr 21 18 0 = .

A1 -Equations for currents in the model
Gating scheme for the 6-state Markovian model for the DHP-sensitive L-type Ca 2+ release channel: The open state O3 dhpr accounts for the increased tail current produced as the result of a large depolarization. Expressions for rate constants: k dhpr where [Ca 2+ ] NaCa , [Na + ] NaCa are Concentrations at the mouth of the NaCaexchanger.
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 .  (85) where [Ca 2+ ] ryr is the Ca 2+ concentration at the mouth of the RyR channel on the dyadic side. In the presence of caffeine (CF, concentration in μM):