Modeling CICR in rat ventricular myocytes: voltage clamp studies
© Krishna et al; licensee BioMed Central Ltd. 2010
Received: 17 June 2010
Accepted: 10 November 2010
Published: 10 November 2010
The past thirty-five years have seen an intense search for the molecular mechanisms underlying calcium-induced calcium-release (CICR) in cardiac myocytes, with voltage clamp (VC) studies being the leading tool employed. Several VC protocols including lowering of extracellular calcium to affect Ca2+ loading of the sarcoplasmic reticulum (SR), and administration of blockers caffeine and thapsigargin have been utilized to probe the phenomena surrounding SR Ca2+ release. Here, we develop a deterministic mathematical model of a rat ventricular myocyte under VC conditions, to better understand mechanisms underlying the response of an isolated cell to calcium perturbation. Motivation for the study was to pinpoint key control variables influencing CICR and examine the role of CICR in the context of a physiological control system regulating cytosolic Ca2+ concentration ([Ca2+] myo ).
The cell model consists of an electrical-equivalent model for the cell membrane and a fluid-compartment model describing the flux of ionic species between the extracellular and several intracellular compartments (cell cytosol, SR and the dyadic coupling unit (DCU), in which resides the mechanistic basis of CICR). The DCU is described as a controller-actuator mechanism, internally stabilized by negative feedback control of the unit's two diametrically-opposed Ca2+ channels (trigger-channel and release-channel). It releases Ca2+ flux into the cyto-plasm and is in turn enclosed within a negative feedback loop involving the SERCA pump, regulating[Ca2+] myo .
Our model reproduces measured VC data published by several laboratories, and generates graded Ca2+ release at high Ca2+ gain in a homeostatically-controlled environment where [Ca2+] myo is precisely regulated. We elucidate the importance of the DCU elements in this process, particularly the role of the ryanodine receptor in controlling SR Ca2+ release, its activation by trigger Ca2+, and its refractory characteristics mediated by the luminal SR Ca2+ sensor. Proper functioning of the DCU, sodium-calcium exchangers and SERCA pump are important in achieving negative feedback control and hence Ca2+ homeostasis.
We examine the role of the above Ca2+ regulating mechanisms in handling various types of induced disturbances in Ca2+ levels by quantifying cellular Ca2+ balance. Our model provides biophysically-based explanations of phenomena associated with CICR generating useful and testable hypotheses.
Contraction of cardiac muscle is triggered by a transient rise in intracellular Ca2+ concentration [Ca2+] myo . Sarcolemmal (SL) membrane depolarization triggers Ca2+ influx from the extracellular medium by opening dihydropyridine (DHP)-sensitive L-type Ca2+ channels. Following diffusion across a small sub-membrane dyadic space, this influx activates ryanodine receptors (RyRs) controlling ryanodine-sensitive Ca2+ release channels in the junctional portion of the sarcoplasmic reticulum (jSR). Fabiato and Fabiato  named the process calcium-induced calcium release (CICR). Ca2+ subsequently diffuses from the dyadic space into the cytosol. Ultimately, intracellular Ca2+ concentration [Ca2+] myo is returned to resting levels by combination of: (a) Ca2+ buffering in the dyadic space and cytosol; (b) sequestration of Ca2+ by sarcoplasmic/endoplasmic reticulum Ca2+-ATPase (SERCA)-type calcium pumps lining the longitudinal portion of the sarcoplasmic reticulum (LSR); and (c) Ca2+ extrusion from the cytosol by Na+/Ca2+ exchangers and Ca2+-ATPase pumps on the sarcolemmal membrane.
CICR in cardiac muscle exhibits both graded behavior and a high gain. Graded behavior refers to the obser-vation that SR Ca2+ release is proportional to the influx of trigger Ca2+ , whereas high gain indicates that the SL trigger current elicits a high SR Ca2+ release flux. Graded Ca2+ release with high gain is somewhat paradoxical according to Stern , in that the positive feedback inherent in such high-gain systems tend to produce regenerative, nearly all-or-none release rather than graded release. Several deterministic models have been developed to explain excitation-contraction (E-C) coupling [4, 5], but none of them can explain the mechanism of graded release at high gain over a wide range of values for sarcolemmal Ca2+ current. Stern  proposed that such a gradation paradox might be explained if the stimulus for Ca2+ release by RyRs were actually the local nanodomains of [Ca2+] generated by nearby L-type channels, rather than the global cytosolic [Ca2+] myo . According to this hypothesis, graded control of macroscopic SR Ca2+ release can be achieved by graded statistical recruitment of individual, autonomous, all-or-none stochastic release events . In these studies, a distributed differential model of high order that included dynamic interactions between large numbers of individual channels was used to demonstrate this concept. However, rather large amounts of computation time are required with distributed stochastic models of this type. Additional models have sought to characterize the Ca2+ release complex, including several [7–9] based on the stochastic release process adopted by Stern et al. These statistical models have solved the graded release problem, however, they too are complicated and computationally very expensive. Other models based on the simplified local control model of CICR developed by Hinch et al.  sought to adopt a lower order description of the E-C coupling process [11, 12] by making an approximation of rapid equilibrium in the dyadic space. The latency from onset of Ca2+ entry via the I Ca,L channel to triggered SR Ca2+ release is known to increase with decrease in the magnitude of I Ca,L,TT , the modeling of which is made possible by considering Ca2+ diffusion in the dyadic medium. These models [11, 12] also approximate the SR as a single volume compartment with no distinction between junctional versus the longitudinal (network) SR compartments. However, recent work  points towards the important role of the Ca2+ refilling rate from the network to junctional SR in controlling RyR release termination via the luminal sensor. Shiferaw et al.  developed a computationally tractable model of Ca2+ cycling to represent the release of calcium from the SR as a sum of spatially localized events that correspond to Ca2+ sparks, assuming the recruitment rate of Ca2+ sparks is directly proportional to the whole-cell I Ca,L current. This assumption overlooks the complex calmodulin mediated interaction (calcium dependent facilitation (CDF) and calcium dependent inactivation (CDI)) of the I Ca,L channel with calcium in its vicinity. It also demands a large amount of computation.
Numerically, distributed as well as statistical models tend to be computationally expensive due to the in-herent repetition involved in the computation. In a spatially distributed model, simultaneous solution for dynamics in identical compartments distributed in space would amount to a large computational cost. In statistical models inference is drawn based on multiple runs of identical events which translate into a prolonged simulation time. Hence, these models are cumbersome to implement, particularly in larger multiple-cell simulations. Consequently, we consider a deterministic approach to the characterization of CICR. Specifically, we develop a lumped model of the Ca2+ release complex that includes: (a) a sub-sarcolemmal dyadic cleft space separating the SL and jSR membranes; (b) a single DHP-sensitive Ca2+ channel on the SL membrane; and (c) a single equivalent Ry-sensitive channel arranged symmetrically on the opposing jSR membrane that represents the output of a local cluster of Ry-sensitive channels facing the DHP-sensitive channel. Based on morphological data compiled by Bers , we further assume that each ventricular cell contains 10,000 of these dyadic Ca2+ release units, and that they are associated with the fraction of the SL membrane that is coupled with the jSR. That is, we partition the sarcolemma into free and dyadically coupled SL membrane, and associate each with a different fluid compartment: the cell cytosolic medium in the case of the free SL membrane, and the dyadic cleft space medium in the case of the dyadic-coupled fraction. In a sense, we build on Stern's  local domain concept by considering the aforementioned local nanodomains identical, but focusing on the nonlinear dynamics of the two different types of Ca2+ channel in the dyadic coupling unit. Our deterministic model although is very descriptive, is computationally tractable and has a run time of 21 sec (including recording of 73 variables of type double on a data file) for 1 cycle of 4 Hz voltage clamp stimulation.
Rat ventricular myocytes were prepared from 200-300 g male Sprague Dawley rats by dissociation with collagenase, as previously described . All experiments were performed under conventional whole cell recording conditions with a List EPC-7 patch clamp, recording fluorescence from nearly the entire cell, as described by Fan and Palade . Recordings from an individual cell were rarely extended beyond 10 min in order to reduce as much as possible both escape of dye from the cell and Ca2+ current rundown. External solution in the bath was normal Tyrode (1 mM Ca2+) with Cs+ substituted for K+ for purposes of blocking inward rectifier K+ currents. The internal solution in the pipette contained Cs aspartate supplemented with 20 mM CsCl, 3 mM Na 2 ATP , 3.5 mM MgCl 2 and 5 mM HEPES. Holding potential used was -40 mV.
All simulations and analysis were performed on a 2.8 GHz Intel® Core™2 Duo CPU-based computer using Microsoft Windows XP operating system. To find the parameters involved in the 6 state Markovian model for I Ca,L , a non-linear least-squares method  was used to obtain the solution of the system of non-linear ordinary differential equations. Specifically, we have employed an algorithm given by Lau . The numerical integration scheme used to solve the full set of forty two 1st-order differential equations describing the dynamic model was the Merson-modified Runge-Kutta 4th-order method [20, 21] with a conservative fixed time step, chosen small enough to allow the local truncation error to be of fourth order. The explicit finite difference scheme was used to numerically solve the Laplacian equations of Ca2+ diffusion in the cleft space. Detailed numerical methods are similar to those presented by Smith et al. . The results were visualized using Matlab by Mathworks and Origin by Microcal Software.
Our objective was to develop a model of the rat ventricular cell which could be used to explain Ca2+ signaling at the nanoscale level of the dyad and integrate the contributions of many dyads to produce a Ca2+ transient and continuous Ca2+ balance at the whole-cell level. Therefore, we start with a broad discussion of the elements of the DCU and its Ca2+ supply (the jSR), and continue with a progressively more detailed description of the whole cell model. It is important to note that all Ca2+ concentrations discussed in the model pertain to unbound Ca2+ unless specified.
Surface area of various plasma membranes in the cell
Surface area of external SL
11.4 × 103 μm2
Surface area of T-tubule
5.52 × 103 μm2
Surface area of total SL (including external SL and T-tubule)
16.9 × 103 μm2 (*)
Surface area of junctional external SL
0.846 × 103 μm2
Surface area of junctional T-tubule
2.54 × 103 μm2
Surface area of total junctional plasma membrane
3.39 × 103 μm2
Surface area of junctional SR
6.99 × 103 μm2
Surface area of longitudinal SR
36.8 × 103 μm2
Surface area of total SR
43.8 × 103 μm2
Parameters used to model sub-cellular morphology
Number of dyadic units
5.3581 × 10-2 nL
Longitudinal SR volume
1.1776 × 10-3 nL
Total junctional SR volume
1.104 × 10-4 nL
Step size in the 'r' direction
Diameter of the cylindrical cleft space in the 'r' direction
Step size in the 'z' direction
Length of the cylindrical cleft space in the 'z' direction
Volume of a unit dyadic space
1.91 × 10-9 nL
Channel and Exchanger Distribution
Parameters used to model ion transport across the sarcolemmal membrane
96485 coul · mol-1
Ideal gas constant
8314 mJ · mol-1 · K-1
[Ca 2+ ] o
Extracellular Ca2+ concentration
[Na + ] o
Extracellular Na+ concentration
[Cs + ] o
Extracellular Cs+ concentration
Z Na , Z Cs
Valence of Na+ and Cs+ ions
Z Ca , Z Ba
Valence of Ca2+ and Ba2+ ions
Permeability of L-Type calcium channel to Ca2+
6.7367 × 10-9 μL · s-1
Permeability of L-Type calcium channel to Na+
8.0355 × 10-11 μL · s-1
Permeability of L-Type calcium channel to Cs+
6.2088 × 10-11 μL · s-1
Dissociation constant for allosteric Ca2+ activation
125 × 10-6 mM
Dissociation constant for extracellular Ca2+
Dissociation constant for intracellular Ca2+
Dissociation constant for extracellular Na+
Dissociation constant for intracellular Na+
Maximum Na+/Ca+ exchange current
Half saturation constant for the SL Ca2+ pump
Maximum sarcolemmal Ca2+ pump current
Dissociation constant for extracellular Cs+
1.5 × 103 μM
Dissociation constant for intracellular Na+
2.14 × 105 μM
Maximum Na+/Cs+ pump current
Maximum background Na+ current conductance
Mean access resistance of the tubular system
The SR Fluid Compartment
Parameters used to model intracellular ion transport
rate of PLB dp phosphorylation
rate of PLB p dephosphorylation
Maximal binding/release rate of Ca2+ from cytosol to SERCA
Maximal binding/release rate of Ca2+ from SERCA to LSR
Total amount of SERCA
Phospholamban to SERCA ratio
PKA ac t
Relative regulatory activity of PKA
Time const. for transfer of Ca2+ from LSR to jSR
7.0 × 10-3 s
Time const. for transfer of Na+ from dyad to myoplasm
1.0 × 10-3 s
Time const. for transfer of Cs+ from dyad to myoplasm
6.0 × 10-3 s
Permeability of the RyR channel
1.714 × 10-7 μL · s-1
Diffusion constant for Ca2+ in the dyadic space
Density of high-affinity Ca2+ binding sites
Density of low-affinity Ca2+ binding sites
Half-saturation value of low-affinity Ca2+ binding site
Half-saturation value of high-affinity Ca2+ binding site
Intracellular Mg2+ concentration
[fluo 3] tot
total concentration of indicator dye
association rate of Ca2+ binding to dye fluo-3
80 μM -1 · s-1
dissociation rate of Ca2+ binding to dye fluo-3
The Dyadic Coupling Unit (DCU)
In describing this functional unit it is necessary to provide progressively more detailed descriptions of the component elements of the individual dyad, particularly the geometrically opposed DHP and Ry-sensitive Ca2+ channels, as well as the geometry and buffering properties of the small dyadic space. As will be shown, we have described the dynamics of both types of Ca2+ channels using Markovian state models (Figure 4) which include features such as Ca2+ mediated channel inactivation, a graded CICR process with a "calcium gain" of approximately 6-7, and two-dimensional Ca2+ diffusion within the dyadic space. Crank  discusses diffusion problems in a two-phase heterogeneous medium and shows that diffusion through a system of barriers (RyR feet structures in the dyadic cleft space) can be approximated by diffusion in the same region without barriers but with a reduced effective diffusion coefficient. We hence take this approach in modeling the Ca2+ diffusion by solving the 2-D Laplacian equation Appendix A3 (Equations 147-150) in the DCU without explicitly accounting for local potential fields. The DHP-sensitive i Ca,L,TT channel brings in trigger Ca2+ (0.1 pA which is of the same order as measured by Wang et. al. ) causing a sparklet (a local increase in Ca2+ concentration at the mouth of the channel). This trigger Ca2+ causes a release from a cluster of opposing RyR channels, causing a spark. This combined release from a cluster of RyR channels causing a spark is represented as the release from a unitary RyR channel (i RyR ) in our model (shown in Figure 6). The characteristics of elemental Ca2+ release from a unitary RyR channel in our model agrees with data in terms of amplitude which is of the order of 3 pA (reported by Cheng et al.  and Blatter et al. ) and duration (full duration at half maximum (FDHM)) which is of the order of 50 ms (reported by Zima et al. ).
The single DCU in our model represents the lumped activity of a large number of individual dyads (e.g. 10,000), and it is charged with the task of forming the cytosolic Ca2+ transient (hence mechanical contraction) each beat of the cardiac muscle cell. In response to tonic application of voltage clamp pulses, the DCU strongly depends on an adequate supply of Ca2+ from the SR. The measurements of Diaz et al.  show that, although trigger current may be supplied regularly by tonic voltage clamp pulses, there is an inherent steady-state dependence of the magnitude of Ca2+ release on the particular value of SR Ca2+ content (i.e., there is a relationship between SR Ca2+ content and peak [Ca2+] myo ; Figure 4 Diaz et al. ). This plot gives us a glimpse of the input-output relationship of the dyadic coupling unit and indicates that SR Ca2+ content is an important controlling variable for the CICR process implemented by the DCU model.
L-Type Ca2+ Current
Rate constants used to model Ca/CaM binding and CaM buffering
3.5 × 10-4‡
1.4 × 10-6‡
Rate constants used to model Ca/CaM/CaMKII interactions
k PP 1
k mPP 1
3.35 × 10-3
2.2 × 10-3
Rate constants used to model Ca/CaM/CaN interactions
1.3 × 10-3
Our previous study  was focused on the characterization of the I Ca,L channel under conditions of low [Ca2+] in the dyadic space and myoplasm. In fact, Ca2+- release from the jSR was blocked by administering a relatively high dose of ryanodine (20 μ M) in all experiments. Therefore, to study the additional influence of SR Ca2+ release on I Ca,L , we modified the original I Ca,L model to better characterize the process of Ca2+-dependent inactivation. In the modified I Ca,L model, majority of the structure for the 6-state dynamic scheme remains the same, however changes have been made to the voltage-dependent inactivation rate function (), and the Ca4CaM dependent inactivation rate function (). The specific formulas for the Markovian state equations and the modified and functions used in this study are given in Appendix A1 (Equations 4-50). Our adjustments consist of reducing the contribution from voltage-dependent inactivation process and strengthening the Ca4CaM-dependent inactivation process. The rate constants and are constructed in order to provide a re-excitation window during the return of the clamp voltage to the resting potential. With these adjustments, the Markovian state description for I Ca,L can provide good fits to measured I Ca,L data under both test conditions (presence and absence of ryanodine -sensitive Ca2+ release) as well as produce tail currents during repolarization from large clamp voltages (≥ 40 mv).
Calcium Buffering in the Dyadic Space
Previous modeling work by Post et al. , and Post and Langer  considered the effect of Ca2+ binding sites on the inner sarcolemmal leaflet. Following these authors, we included low-affinity (k d = 1.1 mM) and high-affinity (k d = 1.3 μ M) Ca2+ binding sites on SL wall boundary of cylindrical dyadic space. The presence of these membrane bound sarcolemmal Ca2+ binding sites in our model has significant physiological implications, in that it prevents local dyadic Ca2+ concentration near the "mouth" of DHP-sensitive I Ca,L channel from becoming excessively high. Allowing such a condition to occur can cause a reversal of the I Ca,L current, which does not occur physically during normal jSR release. Addition of these SL Ca2+ binding sites does not significantly slow the build-up of Ca2+ within the dyadic space, nor is the Ca2+-induced Ca2+ release mechanism affected .
Ca2+ Release Channel
The gating characteristics of the Ry-sensitive release channel are not only modulated by the dyadic Ca2+ concentration at its mouth but also the jSR Ca2+ concentration via the luminal sensor. Several RyR gating schemes have been deduced from isolated RyR currents measured in lipid bilayers, including: (1) the 4-state scheme developed by Keizer and Levine  with 4 Ca2+ binding for activation and 3 Ca2+ binding for inactivation to explain the "adaptation" of the RyR observed by Gyorke and Fill ; (2) the 6-state scheme suggested by Zahradnikova and Zahradnik , which allowed opening the RyR channel upon binding of a single calcium ion; and (3) the Markov model proposed by Keizer and Smith , which can be dynamically switched among the six, five and four-state representations during the simulation as Ca2+ levels vary. Stern et al.  demonstrated that none of these schemes yielded stable local control of SR release, even with extensive adjustment of parameters. Stern et al.  further reported that all schemes resulted in models that manifested local instability, as indicated by failure of release to terminate after activation, or global instability caused by spontaneous activation by resting [Ca2+] myo . Since many of the kinetic gating schemes derived from lipid bi-layer data fail to support stable E-C coupling in simulations, he concluded that the RyR gating process in situ may differ considerably from that in bi-layers.
Our gating scheme is patterned after the release channel used in the model of Stern, Pizarro and Rios  (Figure 4B), where the channel is assumed to have four states: rest (closed), activated (open), inactivated (closed) and refractory (closed). The activation gate is opened by the simultaneous, cooperative binding of two Ca2+ ions, whereas inactivation depends on the binding of a single Ca2+ ion. Four additional features have been built into our model: (a) the crucial role of the luminal SR Ca2+ sensor (Figure 5) in assisting the inactivation of the RyR channel is modeled via the dependence of all the rate constants in the 4-state RyR model, on the degree of interaction between the RyR and the proteins triadin and junctin (Figure 4B); (b) CaMKII act dependent enhancement in RyR release [47–50] is modeled via the rate functions , , and ; (c) a stronger Ca2+ dependent inactivation of the RyR channel (at a fixed depolarization level), is adopted to reflect recent observations that the inactivation of the RyR depends on the high local Ca2+ concentration consequential to their own Ca2+ release ; and (d) inactivation of the RyR at different depolarization levels is made dependent on local Ca2+ concentration (i.e., [Ca2+] ryr at the "mouth" of the RyR channel on the dyadic side) as per Wier et al.  and Zucchi et al. . Our initial studies indicated that the repriming rate () in the RyR gating scheme of Stern et al.  was quite large, which can lead to a saturated open probability of the Ca2+ release channel. This occurs during the later phase of channel inactivation, where the large value of tends to reactivate the channel, resulting in saturated Ca2+ release. Therefore, we utilized a value of that is 10% of that used by Stern et al. , and further assumed that the unitary permeation flux of the jSR release channel is proportional to jSR luminal Ca2+ concentration ([Ca2+] jSR ). The specific equations for the Ca2+ release model are given in Appendix A1 (Equations 73-92).
Luminal RyR sensor
The RyR Ca2+ release is modulated by a multi-molecular Ca2+ signalling complex which is localized to the junctional SR [54–56]. This complex consists of the ryanodine receptor (RyR) which functions as a Ca2+ conducting pore [57, 58], calsequestrin (CS) which acts as the Ca2+ binding protein [59, 60], and the junctional SR transmembrane proteins triadin  and junctin . It was known previously that the proteins triadin and junctin anchor calsequestrin to the ryanodine sensitive receptor . More recently, it has been observed that the protein-protein interactions between triadin, calsequestrin and RyR modulate sarcoplasmic reticulum calcium-release in cardiac myocytes . Triadin and junctin are structurally homologous proteins , but the functional differences in their roles are unclear at present. Therefore we have refrained from modeling these proteins separately with unique roles, and will hereafter only mention triadin. The RyR model proposed by Shannon et al.  takes a heuristic approach towards free SR Ca2+ concentration dependent luminal control of the RyR channel. Our detailed biophysical model of the luminal sensor is strongly based on the recent findings of Terentyev et al.  which uncovers complex [Ca2+] jSR ) dependent, CS mediated mechanistic interaction of the protein triadin with the RyR channel.
The triadin protein facilitates SR Ca2+ release by sensitizing the RyR to activation by the trigger current I Ca,L,TT . This is incorporated in our model by allowing the rate constants in the 4-state Markovian model for the RyR (Figure 4B-ii) channel to be functions of the activated state A 1 ls , which represents the degree of binding between triadin and RyR, in the 6-state model for the luminal sensor (Figure 4B-i).The degree of triadin assisted anchoring of CS to the RyR channel [61, 62] is denoted by the state I 2 ls . Triadin is also known to exist in a form bound to CS alone  denoted by the state I 3 ls and in its free unbound form represented as I 4 ls in the model. CS which is known to exist in a Ca2+ bound form modeled by B 6 ls also modulates SR Ca2+ release by influencing the open probability of the RyR channel [65–67], via interactions with Triadin [54, 66, 25]. The degree of unbound CS is denoted by B 5 ls in the model.
During the diastolic period the relatively large concentration of available free Ca2+ in the SR results in most of the CS being bound to Ca2+, decreasing the degree of interaction between CS and triadin. This enables a strong interaction between the available unbound triadin and the RyR channel increasing its propensity for trigger Ca2+ induced SR release. Our model incorporates this property by facilitating movement of states towards A 1 ls and B 6 ls in the presence of large SR Ca2+ concentration. Following SR Ca2+ release, a reduced luminal Ca2+ concentration in the jSR causes an increase in the amount of free calsequestrin (B 5 ls in the model) available to bind with triadin. This results in a decrease in the extent of interaction between triadin and RyR (A 1 ls in the model), thus inhibiting the RyR channel and leading to robust termination of SR Ca2+ release in cardiac myocytes . This release termination mechanism is incorporated in our combined RyR-Luminal sensor model (Figure 4B and Appendix A1, Equations 73-92).
Parameters used in the luminal sensor model
Ca2+ Buffering in Myoplasm and SR
Ca2+ buffers play an important role in sequestering a fraction of the total Ca2+ released during E-C coupling and contraction. These buffers include: (a) calmodulin (CaM), which is assumed to be uniformly distributed in the myoplasm and dyadic space; (b) troponin in the bulk myoplasm; and (c) calsequestrin (CS) in the jSR. The dyadic space has been shown to be accessible to calmodulin (CaM), but mostly inaccessible to fluorescent dyes [70, 71]. Therefore, we consider the dyadic space filled with calmodulin, but not fluo-3. This provides a more direct calcium communication pathway between the DHP and Ry receptors. The rate constants for Ca2+ binding to calmodulin were based on a model from , whereas those for Ca2+ binding to troponin were taken from Potter and Zott . Rate constants used to describe Ca2+ binding to calsequestrin were based on the study of Cannell and Allen , whereas those for Ca2+ binding to troponin-Mg complex were adopted from Lindblad et al. .
It is well recognized that fluorescent indicator dyes introduced into the cytosol also act as Ca2+ buffers , even at submillimolar concentrations . In our simulations, we have used 100μ M fluo-3. We assume that the quantity observed experimentally as a "calcium concentration" signal is actually the calcium complexed with fluo-3, or [CaF3]. The differential equation describing the change in [CaF3] with time is given in Appendix A2 (Eq. 138) which follows from Shannon et al. .
Cytosolic Ca2+ is pumped into the LSR (Figure 1A) by a Ca2+-ATPase, which lowers [Ca2+] myo and helps to induce relaxation in cardiac muscle. The transport reaction involves two Ca2+ ions and one ATP molecule , and it is represented by the description given for I cyt,serca and I serca,sr in Appendix A1 (Equations 51-63). The model used for the uptake pump is adopted from Koivumaki et al.  and takes into account both the forward flux of Ca2+ from the cytosol to the LSR lumen and the backward flux from the LSR to the cytosol along with the Ca2+ buffering action of the SERCA protein. The phospholamban (PLB) to SERCA ratio has been fixed to 1.0, assuming almost equal availability of both the proteins. CaMKII act affects the SERCA pump via direct phosphorylation assisting in enhancement of SR Ca2+ transport by increasing the pumping rate  and indirectly via phosphorylation of PLB  relieving the inhibition caused by PLB on the SERCA pump in turn increasing the sensitivity of the pump for Ca2+ uptake. These two effects are modeled by allowing the rate constants for Ca2++ binding to/release from the SERCA pump as well as the rate constant for phosphorylation of PLB to be a function of CaMKII act in the cytosol. The activating role of CaN in modulating the SERCA pump [80, 81] is accounted for in our model by allowing the rate constant for phosphorylation of PLB to be dependent on available CaN act in the cytosol. The relative regulatory role of the enzyme protein kinase A(PKA) is fixed at 0.1, as its modulatory effect is beyond the scope of this study. The current I cyt,serca dictates the transport of Ca2+ between the cytosol and the SERCA protein. Similarly, the current I serca,sr dictates the transport of Ca2+ between the SERCA protein and the LSR. The difference in these Ca2+ currents accounts for the Ca2+ buffered by the SERCA protein. The jSR is subsequently refilled by Ca2+ diffusion from the LSR. The differential equations describing the Ca2+ balance and particularly the transfer between two separate SR compartments (jSR and LSR) are provided in Appendix A1 (Eq. 72).
Ca2+-Extrusion via Sarcolemmal Ca2+ Pump
Although the sarcolemmal Ca2+-pump has a high affinity for [Ca2+] myo , its transport rate is far too slow for it to be an important factor in Ca2+ fluxes during the cardiac cycle. It might, however be more important in long-term extrusion of Ca2+ by the cell. Our model of the plasma membrane Ca2+ ATPase pump current is adopted from Sun et al. . We have used a constant value of half activation constant k mpca = 0.5μ M in our model on the basis of measurements by Caroni et al .
Ca2+-Extrusion via Na+/Ca2+ Exchanger
In mammalian cardiac cells, it is generally accepted that the Na+/Ca2+ exchanger has a stoichiometry of 3Na+:1Ca2+ . I NaCa (again, the combination of I NaCa,TT and I NaCa,SL ) is important in removing Ca2+ during twitch relaxation, in competition with I up . A simple thermodynamic Na+/Ca2+ exchanger current model  may be sufficient to predict the direction of Ca2+ transport by Na+/Ca2+ exchange and the driving force, however the amplitude is subject to kinetic limitations (depending on substrate concentrations). A more comprehensive Na+/Ca2+ exchanger current equation  is adopted in our model, which includes factors for allosteric Ca2+ activation and the transport for Na+ and Ca2+ inside and out. The maximal flux through the exchanger V max is estimated to ensure that the Ca2+ ion transport (which is voltage dependent) via the Na+/Ca2+ exchanger matches the influx of Ca2+ via the I Ca,L channel [86, 87], maintaining whole cell Ca2+ homeostasis.
The Na+/K+ pump, helps in maintaining homeostasis of the intracellular Na+ ion concentration. ATPase activity powers the pump, as it generates an outward Na+ flux and an inward K+ flux with a Na+ to K+ stoichiometry of 3 to 2 . However, in our experimental protocol, external solution in the bath was normal Tyrode (1 mM Ca2+) with Cs+ substituted for K+ in order to block the inward rectifier K+ current. The internal solution in the pipette contained Cesium aspartate supplemented with 20 mM CsCl, 3 mM Na2ATP , 3.5 mM MgCl2 and 5 mM HEPES.
Activation of the electrogenic sodium pump in mammalian non-myelinated nerve fibres , skeletal muscle  and rat brain cells  by Cs+ is reported in the literature. It is observed that the late effects of reducing extracellular K+ concentration ([K+] o ) to 0 mM in mammalian cardiac muscle can be prevented by including appropriate concentrations of other activator cations of the Na+/K+ pump such as Cs+ in the 0 mM [K+] o bathing solution . Monovalent cations (including Cs+) were also added to K+ free bathing solution to reactivate the sodium pump in guinea-pig ventricular myocytes . The effectiveness of Cs+ as an external cation in activating the electrogenic sodium pump is known to be lesser than potassium .
We have represented I NaCs in our model by the expression for Na+/K+ pump formulated by Linblad et al.  replacing K+ ion concentrations with the Cs+ ion concentration. While ensuring whole cell Na+ ion balance, the peak Na+/Cs+ pump current is modified to be one-sixth to account for the decreased potency of the cation Cs+ in activating the pump. The voltage-dependence of I NaCs is adopted from the data on Na+/K+ pump from Hansen et al. .
The DCU is a fundamental element in the mechanism of CICR. The sequence of events resulting in CICR is triggered by the Ca2+ entering through the I Ca,L channel. The characteristics of this channel are hence examined in detail to understand its voltage and Ca2+ dependent behavior. The ability of the I Ca,L channel to facilitate graded RyR release is noted. The high gain associated with RyR release is quantified. This is followed by a study of the properties of the RyR channel with a particular emphasis on its inactivation mechanism. The Ca2+ mediated interaction between these channels is investigated in detail. This is followed by an effort to understand the role of [Ca2+] jSR in CICR. In particular, the relationship of the peak cytosolic calcium transient and the SR Ca2+ content is examined. The effect of modulatory agents like caffeine and thapsigargin on SR release/uptake is studied to understand the significance of these mechanisms in facilitating a normal CICR.
L-type Ca2+ current (ICa,L)
Measured data obtained using the same voltage clamp pulse under all three conditions (Control, Ryan-odine application and Barium substitution) is shown for comparison (data obtained from Dr. Palade's lab was pre-processed to eliminate transients produced by changing clamp voltage in order to obtain model fits using a non-linear least-squares method ). The degree of I Ca,L channel opening is known to be enhanced in the presence of CaMKII act  which is known as Ca2+-dependent facilitation (CDF). According to Bers et al.  the positive regulation of the channel by CaMKII act requires Ca2+ influx (it is not seen when Ba2+ is the charge carrier as in Figure 7D and is more strongly apparent when local Ca2+ influx is amplified by SR Ca2+ release as in Figure 7A). Increase in CaMKII act seems to play an important function in enhancing I Ca,L peak amplitude , suggesting a critical role for Ca2+-dependent facilitation. Our model generated results also indicate that a stronger CDF (in Figure 7A) in the presence of RyR Ca2+ release causes an enhancement in I Ca,L peak amplitude (compare the model fits to data in Figure 7A (presence of RyR release) with Figure 7B (no RyR release) and Figure 7C (Ca2+ substituted by Ba2+)). It is important to note that the peak of the I Ca,L current is reached after the RyR open probability reaches its maximum value owing to the faster dynamics of the RyR channel. Our simulations indicate that following the inward current peak, Ca2+-dependent inactivation (CDI) is much faster and dominates the response for the duration of voltage clamp (30 ms <t <80 ms). In addition, a comparison of Figs. 7A and 7B shows that, Ca2+-dependent inactivation (CDI) caused by Ca2+ release from Ry-sensitive channels is much more significant than the self-inhibition produced by Ca2+ influx via I Ca,L itself. The graded behavior of the cytosolic Ca2+ transient is shown in Figure 7C. Voltage dependent inactivation (VDI) is relatively slow compared with CDI and is best seen in Figure 7D where all Ca2+ inactivation effects have been blocked by Ba2+ substitution for Ca2+. Under RyR blockade (Figure 7B), the relatively slow VDI has its major effect during the late phase of the long voltage clamp pulse (e.g., beyond 100 ms) and in a time range where CDI is relatively constant. Quantitative analysis of movement of states (during the depolarizing pulse duration of 50ms) via different pathways in the six state Markovian model shows that 68.24% of the total inactivation of I Ca,L,TT is via the calcium dependent O 2 dhpr -C 4 dhpr pathway, owing to the large Ca2+ concentration which the channel is exposed to in the dyad. In contrast, only 19.21% of the total inactivation of I Ca,L,SL is via the calcium dependent O 2 dhpr -C 4 dhpr pathway, because of the low cytoplasmic Ca2+ concentration in the vicinity of the channel. Considering the total I Ca,L current (I Ca,L,TT + I Ca,L,SL ), around 63.34% of the I Ca,L channel inactivation is via the Ca2+ dependent pathway.
ICa,L,TT - dependent Graded SR-release
The slow rising foot that precedes the rapid upstroke (trace at position 0 in Figure 10B-i and B-ii) is the contribution from a single sparklet (Figure 6), which is the result of Ca2+ release from a single i Ca,L,TT channel in our model (0.1 pA), which is of the same order as reported by Wang et. al.  bringing trigger Ca2+ into a unitary dyad. This L-type Ca2+ channel (LCC) triggers release from a cluster of RyR channels causing a Ca2+ spark (Figure 6), which results in the rapid increase in [Ca2+] dyad that follows the foot (shown in Figure 10B-i and B-ii). As shown in Figure 10B, the latency from the onset of the sparklet foot to the triggered Ca2+ spark increases with decrease in the magnitude of i Ca,L,TT . The SR Ca2+ release flux underlying a typical Ca2+ spark corresponds to approximately 2 pA [27, 28, 13]. From RyR single channel conductance measurements in lipid bilayer studies , a Ca2+ spark translates into a release from around 4 RyR channels (also reported by Blatter et. al. ). This single Ca2+ spark corresponds to a unitary RyR release (2 pA) in our model. Our model does not attempt to reproduce the stochastic kinetics of the single channel LCC-RyR coupling; however it mimics accurately the average behavior of this stochastic process which is the net Ca2+ release flux into the cytosol causing the Ca2+ transient. This on/off stochastic nature of the coupling was used earlier to explain the release termination via stochastic attrition. However, it is now known  that the luminal sensor plays a fundamental role in an active extinguishing mechanism  that effects a robust [Ca2+] jSR - dependent closure of the RyR channel. This mechanism for inactivation of the Ry-sensitive Ca2+ release channel is accounted for in our model.
High Gain of Ca2+ Release
Thus, criterion (1) measures calcium gain observed in the dyadic space, whereas criterion (2) measures gain observed in the cytosolic space.
At a clamp voltage of 10 mv, criterion (1) applied to our simulation yields an integrated flux ratio of 6.39, whereas a gain of 5.75 is calculated using criterion (2). By either method, a CICR amplification factor of approximately 6 is calculated, which is similar to that reported by Stern . This result is also consistent with gain calculations from the measured data of Fan and Palade  on rat ventricular cells. They estimated a gain of approximately 7 using comparisons of the rates of rise of Ca2+ transients in the presence and absence of ryanodine. The model generated results are obtained by using a voltage clamp protocol of a 50 ms step pulse to 10 mv from a holding potential of -40 mv.
The inset in Figure 12A shows the voltage dependence of CICR gain formulated using criterion (1). Our model shows a decline in gain as the clamp voltage is increased from -30 mv to 20 mv [52, 103]. However, any further increase in clamp voltage results in a small increase in gain (Figure 1C, Altamirano et al. ). It is important to note that, with increasing clamp voltage, the decreased ability of the Na+/Ca2+ exchanger (which is co-located  in the dyad) to extrude Ca2+ partially compensates for the declining trigger current, in facilitating SR release and hence assists in increasing CICR gain.
RyR Refractory Characteristics
A dual stimulus protocol (S1-S2) was employed to study RyR refractoriness. The RyR channel was stimulated initially by stimulus trigger current S1, followed after an interval (T2 in Figure 13A) by an identical stimulus current S2. Each of the two stimuli directed into the dyadic coupling unit is the dyadic component (I Ca,L,TT ) of trace 1 (I Ca,L evoked by a 10 mv clamp pulse) in Figure 7A with a peak amplitude of 6.657 pA/pF and a duration (T1 in Figure 13A) of 50 ms. The stimuli S1 & S2 were kept identical to delineate the effects of RyR refractoriness on CICR. The time interval between the two stimuli (T2) was gradually reduced from 250 ms to 50 ms in steps of 25 ms to observe the effects of partial recovery of the RyR channel.
As seen in Figure 13, the decrease in interval between stimuli causes partial recovery of [Ca2+] level in the jSR which manifests in increasing levels of unbound version of the protein calsequestrin, which tends to bind with triadin. As more triadin binds to calsequestrin, there is less left to bind to the luminal side of RyR to activate the channel for CICR. This signaling sequence mediated by the luminal sensor via the interaction of proteins calsequestrin and triadin, allows only a partial SR Ca2+ release based on the degree of SR Ca2+ content recovery. This is an important feature that helps restore the jSR Ca2+ content at the end of every cycle. The RyR channel in our model has an absolute refractory period of 75 ms and a relative refractory period of 250 ms (Figure 13C). These results are consistent with the refractory period measurements by Sobie et al  on rat ventricular myocytes. When the RyR channel is in the absolute refractory period, [Ca2+] RyR drops to a level much below what is caused by a sparklet (trigger Ca2+) and hence robustly avoids re-excitation. However, when the RyR channel is in the relative refractory period, a partial release (Figure 13C) is possible, despite the RyR receptors affinity for release being low.
CICR modulation by the jSR Ca2+ content
From the refractory characteristics of the RyR channel, it is evident that the RyR release is strongly dependent on the jSR Ca2+ content. In fact, the SR Ca2+ release through the RyR channel depends on: (a) the amount of trigger Ca2+ entering by means of the I Ca,L,TT channel (b) the concentration gradient of Ca2+ between the jSR and the mouth of the RyR channel in the dyadic space and (c) the open probability of the RyR channel modulated by the interaction between the luminal sensor and the RyR protein. The SR Ca2+ released also inactivates the trigger current I Ca,L,TT (as shown in Figure 7), hence causing an indirect self-inhibition of RyR release.
Cytosolic peak [Ca2+] dependence on SR Ca2+ content
As seen in traces 1-10 of Figure 16C, during the initial activation of a RyR channel (linear region of phase A in Figure 15), at any specific [Ca2+] ryr the RyR open probability takes a lower value despite a higher steady state diastolic [Ca2+] jSR level. This decreasing slope with increasing diastolic, pre-release [Ca2+] jSR reflects a faster rate of rise in [Ca2+] ryr for a larger pre-release SR content with a very small accompanying increase in the RyR open probability.
Ca2+ Release and its Effect on ICa,L
Adachi-Akahane et al.  investigated I Ca,L -induced Ca2+ release in rat ventricular myocytes in the presence and absence of caffeine, which altered [Ca2+] myo . Using their protocols and data as a guide, we modeled the effect of caffeine on jSR Ca2+ release and its subsequent effect on I Ca,L . Specifically, a Boltzman relationship was adopted to simulate the binding of caffeine to the ryanodine receptor increasing channel conductance, which in turn enhances jSR Ca2+ release. The process is assumed to be dose-dependent, with an open probability of 0.5 at a caffeine concentration of 5 mM (Appendix A1, Eq. 86).
For a given clamp amplitude, a decrease in jSR Ca2+ content decreases Ca2+ release, producing decreased Ca2+ concentration at the dyadic side of the DHP-sensitive channel. Consequently, a smaller amount of calcium-calmodulin complex (Ca4CaM) is activated, leading to a lower degree of Ca2+ dependent inactivation. Thus, I Ca,L peak is enhanced, as shown in Figure 18B, increasing the area under the waveform, corresponding to an increased amount of Ca2+ entering the cell. This inverse relationship between jSR Ca2+ content and the amount of trigger Ca2+ entering the dyadic space for a given clamp voltage is the essence of the feedback mechanism maintaining the efficacy of CICR. This mechanism is operational even in the presence of high concentrations of cytoplasmic Ca2+-buffers, since the small dimensions of the dyadic space allow jSR released Ca2+ to reach its destination (i.e., the DHP receptor) before being captured by the buffers .
Effect of modulation of [Ca2+]o
Calcium balance under conditions of repetitive stimulation
Numerical integration was carried out with respect to steady-state levels, hence only activating currents and leakage currents were considered (background Ca2+ current is neglected). The results of the Ca2+ balance are presented in a manner where relative Ca2+ fluxes can be easily compared. We note that in evaluating the integral of the combined exchanger current I NaCa , the integral is multiplied by a factor of 2 to account for the fact that the exchanger stoichiometry is 3Na+ per Ca2+ (the exchanger transports 1 net charge per Ca2+, whereas all other Ca2+ currents transport 2 charges per Ca2+).
During repetitive long-pulse (800 ms) stimulation, transient Ca2+ fluxes cross the sarcolemmal and SR membranes in both directions (inward (into the cytoplasm) and outward). However, [Ca2+] myo returns to control levels by the end of the long-pulse stimulation protocol. The long-pulse stimulation protocol thus serves as a means of studying the steady-state balance of Ca2+ influx and efflux to and from the cytoplasm, under control conditions. In contrast, this Ca2+ balance is not present in the short-pulse protocol, and [Ca2+] myo is elevated at pulse termination. Our long-pulse simulations show that: (1) the magnitude of individual Ca2+ transients do not change during the 40 pulse sequence; (2) the peak calcium release current (I ryr ) has a constant magnitude; and (3) the occupancy of the calsequestrin Ca2+ buffer in the jSR compartment is reduced by 15% as a result of Ca2+ release during each cycle. Although Ca2+ fluxes cross the sarcolemma and SR membranes in either direction, the integral of each of these currents (indicating charge transfer) over the pulse repetition interval is a constant. Figure 20A shows this well, in that, the charge transfer for each model current (large or small) is a constant, indicating no net loss of Ca2+ nor consecutive-pulse Ca2+ depletion (or augmentation) in the jSR compartment under the long-pulse protocol. In Figure 20, we refer to Ca2+ influx and efflux through the bounding cell membrane (sarcolemma (SL) and T-tubule) as "external fluxes" (i.e., I Ca,L , I NaCa , and I PMCA ), whereas SR membrane Ca2+ fluxes (I ryr , I cyt,serca ) are called "internal fluxes." By convention, inward fluxes are negative (e.g. I Ca,L ) and outward positive (e.g., I NaCa and I PMCA ). Figure 20A shows that the integrated external fluxes sum to zero. The average Ca2+ charge entering the cell via I Ca,L is 603.28 pC (pico-coulomb), and the sum of averaged Ca2+ charges extruded from the cell via I NaCa and I PMCA is 603.28 pC (602.2 pC by I NaCa and 1.08 pC by I PMCA . Figure 20A also shows the magnitudes of the integrated internal fluxes, which sum to zero as well. In calculating inward and outward Ca2+ fluxes relative to the SR lumen, we consider an inward flux (e.g., I cyt,serca ) negative and an outward flux positive (e.g., I ryr ). The integral of I ryr is 751.8 pC, compared with integral of the I Ca,L trigger current (603.28 pC). The calcium gain of CICR calculated in this way is low owing to the increased entry of Ca2+ via the trigger current I Ca,L for a prolonged duration of 800 ms compared to 50 ms previously. The integral of SR Ca2+ uptake is 751.8 pC. These two internal component currents yield flat integral values during long-pulse stimulations and sum to zero.
Decreasing the length of the depolarizing pulse to 100 ms (short pulse) has a pronounced effect on the recovery of contraction (Figure 20B, C and 20D). Compared with the relatively flat amplitude of [Ca2+] myo observed with stimulation by long pulses, short pulse stimulation produces a negative staircase (Figure 20B) in agreement with the data (Negretti et al. ), wherein the peak magnitude decays exponentially to a steady state level of 60% of the initial value. Figure 20B also shows that the model predicts a similar exponential decay in jSR Ca2+ concentration [Ca2+] jSR . Both of these indicators show a net loss of Ca2+ from the cell. Figure 20C and 20D show the transient charge movements associated with the external and internal Ca2+ currents respectively during application of the short-pulse train, which are indicative of an elevated Ca2+ load and resultant Ca2+ imbalance when the pulse train is first applied. Balance is achieved, but at new lower values of [Ca2+] myo and [Ca2+] jSR . For convenient comparison, all internal fluxes are made positive and all external fluxes negative. Figure 20C shows that I NaCa plays a significant role in forming the transient response seen in the short-pulse protocol in extrusion of excess Ca2+ load on the cell. The decreasing amplitude of the Ca2+ transient causes a decreased LSR uptake current (I cyt,serca ), which in turn causes a decrease in [Ca2+] jSR resulting in decreased release (I ryr ). With successive pulses, I cyt,serca and I ryr decline further, mirroring the exponential decline in I NaCa toward balance. At steady state, the external and internal fluxes reach the new equilibrium, where the sum of the integrals of external currents and the sum of the integrals of internal currents are zero.
Long-term calcium stability at higher pacing rates
Ca2+ concentration in the myoplasm
8.1027 × 10-5 mM
Ca2+ concentration in jSR
Ca2+ concentration in LSR
Na+ concentration in myoplasm
Na+ concentration in dyadic space
Cs+ concentration in myoplasm
Cs+ concentration in dyadic space
Fractional occupancy of Calmodulin by Ca2+
Fractional occupancy of Troponin by Ca2+
Fractional occupancy of Troponin by Mg2+ and Ca2+
Fractional occupancy of Troponin by Mg2+
Concentration of Fluo3-Ca2+ complex
C 1 ryr
Closed (resting) state of RyR channel
O 2 ryr
Open (activated) state of RyR channel
1.03668 × 10-9
C 3 ryr
Inactivated state of RyR channel
9.38711 × 10-13
C 1 dhpr
Closed (resting) state of DHPR sensitive Ca2+ channel
O 2 dhpr
Open (activated) state of DHPR sensitive Ca2+ channel
1.499173 × 10-3
O 3 dhpr
Open (activated) state of DHPR sensitive Ca2+ channel
3.300291 × 10-3
C 4 dhpr
Closed (resting) state of DHPR sensitive Ca2+ channel
7.478058 × 10-8
C 6 dhpr
Closed (resting) state of DHPR sensitive Ca2+ channel
7.478058 × 10-8
A 1 ls
Fraction of Tr/J bound to RyR
I 2 ls
Fraction of Tr/J bound to RyR and Calsequestrin
I 3 ls
Fraction of Tr/J bound only to Calsequestrin
B 6 ls
Fraction Calsequestrin bound to Ca2+
Ca bound to the serca protein
36.0637 × 10-5μM
S 2 dhpr
Fraction of IQ Motif bound to Ca4CaM
6.816729 × 10-2
Fraction of unphosphorylated Phospholamban
7.684160 × 10-2
Ca 2 CaM
2 Ca2+ ions bound to C-terminus of CaM
Ca 4 CaM
4 Ca2+ ions bound to C & N terminus of CaM
8.635052 × 10-2
7.563836 × 10-2
Ca 2 CaMB
Ca 4 CaMB
Buffered Ca 4CaM
1.288455 × 10-6
Ca 4 CaN
CaN bound to 4 Ca 2+ ions
2.606246 × 10-4μM
CaN bound to CaM
4.348535 × 10-3μM
Ca 2 CaMCaN
CaN bound to 2 Ca2+ ions and CaM
1.419613 × 10-1μM
Ca 4 CaMCaN
CaN bound to 4 Ca2+ ions and CaM
Ca2+ concentration in the dyadic space
9.012 × 10-5 mM
Fraction of inactive dephosphorylated CaMKII in Ca2CaM bound state
5.527608 × 10-1
Fraction of active dephosphorylated CaMKII in Ca4CaM bound state
3.661260 × 10-1
Fraction of active Thr287-autophosphorylated but CaM autonomous CaMKII
1.314410 × 10-3
Fraction of active Thr287-autophosphorylated but Ca2CaM bound CaMKII
6.277911 × 10-7
Fraction of active Thr287-autophosphorylated but Ca4CaM trapped CaMKII
9.121920 × 10-8
Secondary [Ca2+]myo transients induced by "tail currents"
It is well-established that mammalian cardiac excitation-contraction coupling is mediated by calcium-induced calcium release (CICR). We have developed a comprehensive mechanistic model of CICR under voltage clamp conditions in the rat ventricular myocyte, which includes electrical equivalent circuit models for both the free sarcolemma and that portion involving junctional transmission, as well as, fluid compartment models for several fluid media within the cell (dyadic cleft space, longitudinal sarcoplasmic reticulum (LSR or Ca2+ uptake compartment), junctional sarcoplasmic reticulum (jSR or release compartment), and the cytosolic fluid compartment). An external bathing medium completes our fluid compartment description of the cell Figure 1. The multiple component model is referred to as the "whole-cell model" (Figures 1 and 2). We have probed the mechanisms regulating Ca2+ in the cell. In particular, we have focused on the dyadic mechanisms effecting CICR.
The dyadic controller is a finely tuned coupling device consisting of two opposed Ca2+ channels separated by a small dyadic space (Figure 1B): the sarcolemmal DHP-sensitive "trigger" channel and the Ry-sensitive jSR Ca2+ release channel. The trigger channel is voltage-activated and is driven in our simulations by a voltage clamp pulse which opens the channels, admitting Ca2+ influx (trigger current) to the dyadic space. After diffusion within the dyadic space, trigger Ca2+ affects the Ca2+ dependent open probability of the jSR Ca2+ release channel over a small range of Ca2+ concentrations. Clamp voltage magnitude strongly affects the I Ca,L current and hence the amount of Ca2+ delivered to the dyadic space, but the range over which the Ca2+ concentration at the mouth of the Ca2+ release channel changes (in the absence of RyR release) is quite small (23.42μ M to 8.82μ M for a change in VC from 10 mv to 40 mv). The RyR-mediated Ca2+ release is activated by trigger current, but the release itself is affected by the concentration gradient between the jSR and the dyadic space, and the temporal and refractory properties of the ryanodine receptor (Figure 13). The four-state RyR model is informed regarding supply Ca2+ content by the jSR luminal sensor, a novel feature in our model which characterizes the important protein-protein interactions between calsequesterin, triadin/junctin and the RyR receptor. Triadin/junctin strongly regulates the sensitivity of RyR to trigger Ca2+. Hence the luminal sensor, which is a key element responsible for robust post-release RyR inactivation and refractoriness of the Ry-sensitive Ca2+ release channel, is critical in providing realistic fits to cytosolic Ca2+ transients and an adequate refilling time for the SR Ca2+ stores. Ca2+ release is thus graded with Ca2+ concentration at the mouth of Ca2+ channel in a very sensitive manner with a gain of approximately 7, as is shown in Figure 12.
The sarcolemmal portion of the dyadic membrane defining the dyadic space contains voltage-sensitive Ca2+ channels and deals with changes in the external environment of the ventricular cell (e.g., membrane response to changes in transmembrane potential and chemical signaling agents). It integrates these various stimuli and delivers a trigger current to the small dyadic space. In contrast, the jSR membrane lining the opposite boundary of the dyadic space is concerned with the adequacy of Ca2+ release in CICR. It contains Ry-sensitive Ca2+ channels that require a Ca2+ concentration gradient directed across the channel and into the dyadic space for their operation. Thus, jSR Ca2+ concentration must be maintained within an acceptable range so that calcium is always available for ready release. The relationship between jSR Ca2+ content and peak [Ca2+] myo is shown in Figure 16A. The LSR compartment connects the jSR compartment with the cytosolic compartment. It feeds make-up calcium to the jSR, by using a SERCA pump to actively draw in Ca2+ from the myoplasm. Pumping rate is controlled by the cytosolic Ca2+ concentration [Ca2+] myo as well as the LSR Ca2+ concentration [Ca2+] LSR as part of a mechanism for replenishing and maintaining jSR Ca2+ stores. Our model also incorporates the Ca2+ induced CaM mediated effects of CaMKII and CaN on targets such as the DHP-sensitive I Ca,L channel, Ry-sensitive Ca2+ release channel, as well as the SR Ca2+ ATPase pump. Provision for the effects of phospholamban on SERCA has also been included.
The voltage-gated and chemically gated channels of the dyad are tightly coupled by feedback mechanisms that involve Ca2+ signaling. Although physically separated, the voltage-sensitive Ca2+ channel (Figure 9) as well as the Ry-sensitive release channel (Phase C in Figure 15) are inhibited by increased Ca2+ levels in the dyadic space . In earlier models which lacked the luminal sensor, RyR self-inhibition by dyadic Ca2+ was the only mechanism besides stochastic attrition (both of which are inadequate in effecting a robust RyR channel closure) that was given the role of initiating RyR recovery. This has been tested in our model by artificially clamping all the variables (including concentration at the mouth of the RyR channel) to the levels reached at T2 in Figure 15 despite which the RyR open probability, begins to decline in phase C of Figure 15 showing the self-inhibitory role of high dyadic Ca2+ concentration. Our model predicts (data not shown) that this Ca2+ signaling continues even in the presence of high concentrations of Ca2+ buffering agents in the cytosol (in agreement with the data of Adachi-Akahane et al.  and Diaz et al. ). Tight coupling between the DHP and Ry-sensitive Ca2+ channels within the dyadic space thus preserves the mechanism of CICR under these extreme conditions.
Regulation of cytosolic Ca2+ concentration [Ca2+] myo is evident at the cellular level. It is affected strongly by the voltage and [Ca2+] myo -dependent properties of the sarcolemmal currents I NaCa and I Ca,L , as well as the [Ca2+] myo -dependent properties of the I PMCA and SERCA pumps. We demonstrate a model-generated whole-cell Ca2+ balance, which shows the importance of the Na+/Ca+ exchanger in extruding the Ca2+ that has entered the cell under normal activity, and also any excess that might occur when cytosolic Ca2+ levels rise. The variety of experiments emulated in this study demonstrates quantification of Ca2+ balances for all external and internal Ca2+ fluxes and shows that the model has long-term stability in regulating cytosolic Ca2+, as shown in the 120-sec duration experiments of Negretti et al. (Figure 20) at a pulse repetition rate of 0.33 Hz., and the faster paced stimulation at 4 Hz shown in Figure 14A.
We have examined the dyadic component of the I Ca,L with regard to its Ca2+ -dependent inactivation as a function of jSR Ca2+ content (Figure 17C). Trigger current is voltage-dependent and therefore parameterized by clamp pulse amplitude (Figure 7A). In addition, its inactivation is Ca2+ dependent, especially during Ca2+ release. At a constant depolarizing voltage pulse, the pre-release jSR Ca2+ content dictates the peak of the cytosolic Ca2+ transient (Figure 16A), by controlling peak Ca2+ concentration in individual dyads. As jSR Ca2+ content increases, so does the peak dyadic Ca2+ concentration, and the amount of sarcolemmal Ca2+ influx declines due to greater Ca2+ dependent inactivation. This autoregulatory feedback mechanism helps to establish a stable operating point for jSR Ca2+ content. Transient increases in jSR Ca2+ content bring about increases in Ca2+ release, but reflexly, decrease sarcolemmal Ca2+ influx via I Ca,L . Coupled with other dyadic and extra-dyadic mechanisms (I NaCa ) that decrease [Ca2+] myo and hence Ca2+-uptake via the SERCA pump, jSR Ca2+ content decreases. For decreases in jSR Ca2+ content, the opposite occurs, so that regardless of the sign of the perturbation in jSR Ca2+ content, it tends to stay constant in the steady-state. This is an important feature of the dyadic mechanism that preserves the integrity of CICR. The gain of this feedback system about the operating point is visualized as the slope of the peak Ca2+ transient vs SR Ca2+ content characteristic (Figure 16A) for a given voltage pulse level. Observing this figure, we note that there is a linear operating range beyond which the system gain increases dramatically. Our model results also show that a decrease of jSR Ca2+ content (simulated either by a decrease of Ca2+ uptake into the LSR by thapsigargin or an increase of Ca2+ leak out of the jSR by caffeine) decreases systolic [Ca2+] myo , and hence the model might serve as a useful adjunct in a study of heart failure, where decreased contractility as a result of diminished Ca2+ transients are commonly observed .
This model of a rat ventricular myocyte is limited to Ca2+ related channel, exchanger and pumps (I Ca,L , I NaCa , I PMCA and SERCA pump), while lacking exclusive Na+ or K+ related channels and transporters and is based on data at positive potentials in the range 10 mV ≤ V ≤ 40 mV. It is aimed at mimicking voltage clamp conditions where channels other than calcium are blocked, and it cannot be used to study experiments involving the generation of action potentials. However, its focus on the Ca2+ dynamics allows one to comprehend more clearly the important role of Ca2+ signalling pathways and feedback control systems in maintaining whole cell homeostasis over a prolonged period of time.
A single dyadic space in our model has one representative, lumped DHP-sensitive and Ry-sensitive Ca2+ channel, on opposing sarcolemmal and SR surfaces respectively. This simplified configuration does not conform to detailed structural information regarding the geometrical relationships between DHP and RyR-sensitive Ca2+ channels [115, 55, 116] and thus cannot be used to draw conclusions about this part of the EC coupling process. However, our lumped abstraction forms a functional model of the dyadic coupling unit that can produce accurate predictions of cytosolic Ca2+ transients. The effectiveness of this model is further demonstrated by its ability to accurately characterize the interaction between the DHP and Ry-sensitive Ca2+ channels, including pulse duration dependent termination of release (Figure 11), Ca2+ dependent inactivation of I Ca,L (Figures 7, 9), as well as the wide variety of whole-cell voltage clamp protocols (Figures 17, 18, 19 and 20).
Although our model provides secondary Ca2+ tail transients elicited by I Ca,L "tail currents" (Figure 21), this aspect of the model has not been verified extensively due to paucity of measured VC data showing tail transients over the high voltage range (40 ≤ V ≤ 60) in rat ventricular myocytes. Our I Ca,L tail currents are in general agreement with model generated data at a clamp voltage of 50 mV reported by Greenstein et al. (; Figure 8). The restitution time for the RyR channel in rat ventricular myocytes is believed to be at least 25 ms  and as fast as 150 ms  indicating the wide range of values reported in the literature. We demonstrate model-generated RyR refractoriness by changing the duration of the simulated VC pulse and obtain RyR recovery characteristics that are consistent with measured data for Ca2+ spark restitution in rat ventricular myocytes reported by Sobie et al. . However, our results on this important phenomena, particularly the onset of tail transients with increasing inter-stimuli interval are preliminary and further modeling investigations would benefit considerably from the availability of additional measured data.
We have developed a mathematical model of Ca2+ dynamics under voltage clamp conditions in the rat ventricular myocyte, which is based solidly on experimental data and includes the most extensive description available of a novel feature, namely the luminal Ca2+ sensor in the junctional SR which models the protein-protein interaction between triadin/junctin, calsequestrin and the RyR channel. The luminal sensor imparts the much needed refractoriness to the Ry-sensitive Ca2+ release channel. This element is critical in providing realistic fits to cytosolic Ca2+ transients and an adequate refilling time for the SR Ca2+ stores. Our voltage-clamp simulations demonstrate graded Ca2+ transients with sufficient gain, as well as quantification of Ca2+ balances for all external and internal Ca2+ fluxes. Our model of the dyadic coupling unit (DCU) provides mechanistic explanations of the major input-output relationship for CICR (Figure 16), as well as its modulation by trigger current (clamp voltage). The variety of experiments emulated in this study demonstrates that the model has long-term stability in regulating cytosolic Ca2+, as shown in the 120-sec duration experiments of Negretti et al. (Figure 20) at a pulse repetition rate of 0.33 Hz., and the faster (physiological) paced stimulation at 4 Hz shown in Figure 14A. It also provides biophysically based insights into the molecular mechanisms underlying whole-cell responses to the wide variety of testing approaches used in voltage clamp studies of myocytes that have appeared in the literature over the past two decades (Figures 17, 18, 19 and 20). Thus, the model serves as a platform for the predictive modeling of VC investigations in a number of areas. These include new hypotheses with regards to the under-expression of triadin/junctin resulting in a malfunctioning luminal sensor, which could affect long-term calcium stability of the cell (Figure 14B), and/or changes in the refractoriness of the RyR Ca2+ channel (Figures 13 and 21) affecting the integrity of CICR under a variety of conditions. These are fundamental issues that would benefit from a better mechanistic understanding of deranged calcium signalling in the rat ventricluar myocyte. This study is aimed at providing an initial step towards this goal.
Below is the complete set of equations used in the model.
A1 - Equations for currents in the model
L-Type Ca 2+ current
where k Cs = 0.5 and [Ca2+] dhp , [Na+] dhp , [Cs+] dhp are Concentrations at the mouth of the DHP-sensitive Ca2+ channel.
The corresponding unitary currents i Ca , L , i Na , i Cs are obtained by dividing the above net channel currents by the number of dyadic units, N dyad .
Uptake of Ca2+ from the cytosol into the LSR
Ca2+ pump in SL
where [Ca2+] NaCa , [Na+] NaCa are Concentrations at the mouth of the NaCa-exchanger.
The corresponding unitary current i NaCa is obtained by dividing the above net channel current I NaCa by the number of dyadic units N dyad .
The corresponding unitary current i NaCs is obtained by dividing the above net channel current I NaCs by the number of dyadic units N dyad .
Background Na+ current
Ca2+ transfer from LSR to a single jSR
Ca2+ release from a unit jSR into a single DCU
A 1 ls Fractional occupancy of RyR by Triadin/Junctin
I 2 ls Fractional occupancy of RyR by Triadin/Junctin and Calsequestrin
I 3 ls Fractional occupancy of Triadin/Junctin by Calsequestrin
I 4 ls Free Triadin/Junctin, B5ls - Free Calsequestrin
B 6 ls Fractional Occupancy of Calsequestrin by Calcium
Ca2+ dependent CaM mediated activation of CaMKII and CaN
A2 - Differential equations for buffers used in the model
Intracellular Ca2+ buffering:
A3 - Differential equations for ion concentrations used in the model
Intracellular Ca 2+ concentration: