Open Access

Modeling CICR in rat ventricular myocytes: voltage clamp studies

  • Abhilash Krishna1,
  • Liang Sun2,
  • Miguel Valderrábano3,
  • Philip T Palade4 and
  • John W ClarkJr1Email author
Theoretical Biology and Medical Modelling20107:43

DOI: 10.1186/1742-4682-7-43

Received: 17 June 2010

Accepted: 10 November 2010

Published: 10 November 2010

Abstract

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 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 [1] 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+ [2], 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 [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 Ca2+ current. Stern [3] 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 [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 Ca2+ release complex, including several [79] 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 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 [13], 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 [14] 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. [15] 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 [16], 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 [3] 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.

Methods

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 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.

Computational Aspects

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 [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 Ca2+ 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 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.

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].
Figure 1

Cellular fluid compartments. Figure 1: Cellular fluid compartments. (a) Model configuration showing dyadic space, jSR, LSR, cytoplasm and SL; (b) Inset provides a more detailed description of the dyadic space showing the coupling of the two types of Ca2+ channels (trigger and Ca2+ release channels) via the dyadic fluid medium. Only one representative dyadic coupling unit is shown, however the whole model contains such 10,000 identical units lumped together.

Figure 2

Electrical equivalent circuit for the plasma membrane of a rat ventricular cell. Figure 2: Cm,TT: membrane capacitance of the junctional SL membrane coupled with the dyadic space; Cm,SL: membrane capacitance of the uncoupled free SL membrane; Currents through the uncoupled free SL membrane are (a) ICa,L,SL: L-type calcium current, (b) INa,SL: sodium current through the DHPR channel, (c) ICs,SL: cesium current through the DHPR channel, (d) INaC a,SL: sodium-calcium exchanger current, (e) INaCs,SL: sodium-cesium exchanger current, (f) I PMCA : calcium pump current, (g) INa,b: background sodium current; Currents through the junctional SL membrane coupled with a dyadic space are (h) ICa,L,TT: L-type calcium current; (i) INa,TT: sodium current through the DHPR channel, (j) ICs,TT: cesium current through the DHPR channel, (k) INaCa,TT: sodium-calcium exchanger current, (l) INaCs,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.

Figure 3

Fluid compartment model. Figure 3: Representative fluid compartment model, showing membrane surface area separating different compartments. M FreeSL : free SL membrane; M JunctionalSL : junctional SL membrane; M jSR : junctional SR membrane; M LSR : Longitudinal SR membrane.

Table 1

Surface area of various plasma membranes in the cell

Variable

Description

Value

A Ext.SL

Surface area of external SL

11.4 × 103 μm2

A TT

Surface area of T-tubule

5.52 × 103 μm2

A TotSL

Surface area of total SL (including external SL and T-tubule)

16.9 × 103 μm2 (*)

A JunctExt.SL

Surface area of junctional external SL

0.846 × 103 μm2

A JunctTT

Surface area of junctional T-tubule

2.54 × 103 μm2

A TotJunct

Surface area of total junctional plasma membrane

3.39 × 103 μm2

A JunctSR

Surface area of junctional SR

6.99 × 103 μm2

A LongSR

Surface area of longitudinal SR

36.8 × 103 μm2

A TotSR

Surface area of total SR

43.8 × 103 μm2

*Electrical capacitance of the cell membrane = 169 pF (using 1 μ F/cm2)

Table 1: Surface area of various plasma membranes in the cell. All the above parameters are derived from Bers [16].

Table 2

Parameters used to model sub-cellular morphology

Parameter

Definition

Value

References

N dyad

Number of dyadic units

10000

[8]

V myo

Myoplasmic volume

5.3581 × 10-2 nL

[16]

V LSR

Longitudinal SR volume

1.1776 × 10-3 nL

[16]*

N d y a d V j S R

Total junctional SR volume

1.104 × 10-4 nL

[16]*

Δr

Step size in the 'r' direction

10 nm

Numerical solution

d

Diameter of the cylindrical cleft space in the 'r' direction

400 nm

[119, 8, 120, 121, 6]

Δz

Step size in the 'z' direction

0.76 nm

Numerical solution

h

Length of the cylindrical cleft space in the 'z' direction

15.2 nm

[119, 8, 120, 121, 6]

V cleft

Volume of a unit dyadic space

1.91 × 10-9 nL

--

Table 2: Parameters used to model sub-cellular morphology. Adopted (*), derived (†) or estimated (‡) from the cited sources.

Channel and Exchanger Distribution

Recent research has also shown that besides L-type Ca2+ channels, Na+/Ca+ exchanger activity is also found predominantly in the T-tubules of rat ventricular myocytes [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 Ca2+ into a unitary dyad), i NaCa,TT and i NaCs,TT respectively (Figure 2). The free sarcolemmal component of these same channels are denoted as I Ca,L,SL , I NaCa,SL and I NaCs,SL respectively. I Ca,L,TT is the total current entering through L-type Ca2+ channels via all the dyadic units (N dyad ). We define the total L-type current I Ca,L as the combination of I Ca,L,TT and I Ca,L,SL (i.e., I Ca,L = I Ca,L,TT + I Ca,L,SL ). I Ca,L in our model is mostly (90%) from the L-type Ca2+ current in the T-tubules, since Kawai et al. [24], found L-type current to be highly concentrated (9-fold) in the T-tubules (I Ca,L,TT ) vs. the cell surface sarcolemma (I Ca,L,SL ) of rat ventricular myocytes. We described the I Ca,L channel using a 6-state Markovian model as shown in Figure 4A. The distribution of I NaCa and I NaCs correspond to that of I Ca,L channel in order to ensure Ca2+ and Na+ ion balance.
Figure 4

Calcium channel dynamics. Figure 4: Calcium channel dynamics. (a) Markovian model describing the DHP-sensitive Ca2+ channel, and (b) Markovian model of the Ry-sensitive Ca2+ channel and the luminal SR Ca2+ sensor. Input from the luminal SR Ca2+ sensor modulates the rate constants in the model of Ry-sensitive channel exercising the indirect bias of luminal [Ca2+] jSR on the RyR receptor.

Since our study is focused on voltage clamp testing of Ca2+ transients in rat ventricular cells, we assume that the majority of Na+ and K+ channels are blocked by either the holding potential used (-40 mV) or appropriate blocking agents. Thus, these channels are not modeled, and we assume that only the dihy-dropyridine (DHP) -sensitive Ca2+ channels, the electrogenic pumps, Na+/Ca+ exchangers and Na+/Cs+ pumps expressed in the free and/or coupled SL membranes contribute to the voltage clamp response. Table 3 provides values for various parameters used to model the ion transport across the sarcolemmal membrane.
Table 3

Parameters used to model ion transport across the sarcolemmal membrane

Parameter

Definition

Value

References

F

Faraday's constant

96485 coul · mol-1

--

R

Ideal gas constant

8314 mJ · mol-1 · K-1

--

T

Absolute temperature

290 K

Measured

[Ca 2+ ] o

Extracellular Ca2+ concentration

1.0 mM

Measured

[Na + ] o

Extracellular Na+ concentration

140.0 mM

Measured

[Cs + ] o

Extracellular Cs+ concentration

3.0 mM

Measured

Z Na , Z Cs

Valence of Na+ and Cs+ ions

1.0

--

Z Ca , Z Ba

Valence of Ca2+ and Ba2+ ions

2.0

--

P Ca

Permeability of L-Type calcium channel to Ca2+

6.7367 × 10-9 μL · s-1

[30]*

P Na

Permeability of L-Type calcium channel to Na+

8.0355 × 10-11 μL · s-1

[30]*

P Cs

Permeability of L-Type calcium channel to Cs+

6.2088 × 10-11 μL · s-1

[30]*

K mAllo

Dissociation constant for allosteric Ca2+ activation

125 × 10-6 mM

[85]*

K mCao

Dissociation constant for extracellular Ca2+

1.14 mM

[85]*

K mCai

Dissociation constant for intracellular Ca2+

0.0036 mM

[122]*

K mNao

Dissociation constant for extracellular Na+

87.5 mM

[83]*

K mNai

Dissociation constant for intracellular Na+

12.3 mM

[85]*

V max

Maximum Na+/Ca+ exchange current

776.2392 pA

[85]*

k mpca

Half saturation constant for the SL Ca2+ pump

0.5 μM

[30]*

I ¯ P M C A

Maximum sarcolemmal Ca2+ pump current

1.15 pA

[30]*

K mcs

Dissociation constant for extracellular Cs+

1.5 × 103 μM

[94, 92, 88, 74]

K mna

Dissociation constant for intracellular Na+

2.14 × 105 μM

[94, 92, 88, 74]

I ¯ N a C s

Maximum Na+/Cs+ pump current

147.3 pA

[94, 92, 88, 74]

G Nab

Maximum background Na+ current conductance

0.00141 nS

[30]*

R a

Mean access resistance of the tubular system

20.0

[123]*

Table 3: Parameters used to model the transmembrane currents ICa,L, I NaCa , I PMCA , I NaCs and the background sodium current I Na,b . Adopted (*) or estimated (‡) from the cited sources.

The SR Fluid Compartment

The SR is an intracellular organelle that consists of two lumped fluid compartments (the jSR and LSR) that communicate (Figure 1A; Figure 3). Like the sarcolemma, the bounding membranes of the jSR and LSR are differentiated regarding their ionic current content and degree of coupling with the sarcolemma. With regard to ionic currents, the LSR membrane has a thapsigargin-sensitive SERCA pump for pumping Ca2+ into the LSR lumen against a concentration gradient. In contrast, the jSR membrane contains an outwardly directed ryanodine (Ry)-sensitive channel for Ca2+ release from the jSR to the dyadic space. The jSR fluid compartment contains the Ca2+ binding protein calsequestin as well as the proteins triadin and junctin, which interact with the ryanodine receptor (RyR) and calsequestrin. This co-located configuration of the RyR receptor, along with the proteins calsequestrin, triadin and junctin which exist on the luminal side of the jSR membrane, constitutes a jSR Ca2+ release regulating mechanism called the luminal sensor (Figure 4B; Figure 5). The protein-protein interaction between them plays an important role in regulating the open-state of the RyR Ca2+ release channel [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 Ca2+ sensor. Figure 4B shows a functional diagram of the luminal sensor and its output state is shown connected to the four-state RyR model. Specifically, the sensor adjusts Ca2+ dependent rate functions within the ryanodine receptor model, which affects the open probability P o of the SR Ca2+ release channel.
Figure 5

Schematic of the CICR subsystem. Figure 5: Schematic of the CICR subsystem. The dyadic coupling unit comprises of a DHP-sensitive Ca2+ channel opposing a Ry-sensitive Ca2+ channel. The transmembrane proteins triadin and junction along with calsequestrin mediate the interaction between the luminal Ca2+ and the RyR thus regulating the release of Ca2+ flux from the jSR into the dyadic space.

With regard to coupling, the DHP and Ry-sensitive Ca2+ channels are assumed to be located on opposite sides of the small dyadic fluid space (nanodomain) as shown in Figure 6, and coupled functionally by a CICR mechanism. The dyadic space is assumed to be in fluid communication with the cell cytoplasm via a restricted diffusion region. In contrast, the LSR is not functionally coupled to the sarcolemma, but rather is in contact with the cytoplasm via the SERCA pump (as shown in Figure 3). Table 4 provides values for parameters used to model the intracellular ion transport.
Figure 6

Schematic of the Dyadic nanodomain. Figure 6: Schematic of the LCC-RyR coupling. The dyadic cleft allows Ca2+ signaling between a DHP-sensitive Ca2+ channel opposing a cluster of Ry-sensitive Ca2+ channels (represented as the release channel i RyR in our model). The trigger sparklet due to Ca2+ entering via DHP-sensitive iCa,L,TTchannel causes a spark as a result of release from a cluster of opposing RyR channels. The large amount of Ca2+ released subsequently diffuses into the cytosol.

Table 4

Parameters used to model intracellular ion transport

Parameter

Definition

Value

References

k 12 P L B

rate of PLB dp phosphorylation

6800 s-1

[77]*

k 21 P L B

rate of PLB p dephosphorylation

1000 s-1

[77]*

K cyt,serca

Maximal binding/release rate of Ca2+ from cytosol to SERCA

6250 s-1

[77]*

K serca,sr

Maximal binding/release rate of Ca2+ from SERCA to LSR

6.25 s-1

[77]*

SERCA tot

Total amount of SERCA

47.0 μM

[77]*

PSR

Phospholamban to SERCA ratio

1.0

[77]*

PKA ac t

Relative regulatory activity of PKA

0.1

[77]

τ tr

Time const. for transfer of Ca2+ from LSR to jSR

7.0 × 10-3 s

[30]

τ Na

Time const. for transfer of Na+ from dyad to myoplasm

1.0 × 10-3 s

[30]

τ Cs

Time const. for transfer of Cs+ from dyad to myoplasm

6.0 × 10-3 s

[30]

P ryr

Permeability of the RyR channel

1.714 × 10-7 μL · s-1

[27, 28, 13]

D Ca

Diffusion constant for Ca2+ in the dyadic space

100.0 μm2s-1

[119]*

N h

Density of high-affinity Ca2+ binding sites

200.0 μM

[119]

N l

Density of low-affinity Ca2+ binding sites

16.0 μM

[119]

K l

Half-saturation value of low-affinity Ca2+ binding site

1100.0 μM

[119]*

K h

Half-saturation value of high-affinity Ca2+ binding site

13.0 μM

[119]*

[Mg2+] myo

Intracellular Mg2+ concentration

634 μM

[124]*

[fluo 3] tot

total concentration of indicator dye

100.0 μM

Measured

k f l u o 3 +

association rate of Ca2+ binding to dye fluo-3

80 μM -1 · s-1

[125]*

k f l u o 3

dissociation rate of Ca2+ binding to dye fluo-3

90 s-1

[125]*

Table 4 Parameters used to model the intracellular ion transport I cyt,serca , I serca.sr , I tr , I RyR , and Ca2+ diffusion/buffering. Adopted (*), derived (†) or estimated (‡) from the cited sources.

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 [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 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. [13]) 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. [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 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. [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 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. [29]). This plot gives us a glimpse of the input-output relationship of the dyadic coupling unit and indicates that SR Ca2+ content is an important controlling variable for the CICR process implemented by the DCU model.

L-Type Ca2+ Current

A multiple state characterization of I Ca,L in rat ventricular myocytes has been reported previously by our group [30]. The gating scheme used in this I Ca,L model has an additional high voltage state C 6 dhpr as shown in Figure 4A which is introduced to reproduce I Ca,L tail currents. Upon voltage-dependent activation, the channel achieves the primary open state O 2 dhpr . The degree of opening is enhanced in the presence of activated calmodulin-dependent kinase (CaMKII act ) [3133] which is known as Ca2+-dependent facilitation (CDF). CaMKII has been shown to tether to the I Ca,L channel [34] functioning as a Ca2+ signalling sensor for facilitation. The open probability of the I Ca,L channel is also increased in the presence of activated calcineurin (CaN act ) [35]. These two effects are modeled as shown in Figure 4A, where k 12 d h p r is a function of both CaMKII act and CaN act besides voltage. The gating scheme features two different pathways for inactivation of the open state. The pathway O 2 dhpr C 5 dhpr accounts for voltage-dependent inactivation, whereas the pathway O 2 dhpr C 4 dhpr C 5 dhpr with Ca2+-calmodulin (Ca 4 CaM) dependent rate constant k 24 d h p r accounts for the fast and slow phases of Ca2+-dependent inactivation (CDI) [36]. Transitions to the closed state via both the inactivation pathways are suppressed in the presence of CaMKII act and CaN act to allow CDF. A 2-state Markovian model allows Ca2+ mediated interaction between calmodulin and the I Ca,L channel. As shown in Figure 4A, state S 2 dhpr , which denotes Ca 4 CaM bound to the IQ-motif of the I Ca,L channel, modulates calmodulin dependent Ca2+-induced inactivation. This is in agreement with the findings of Nikolai M. Soldatov [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 and 7 provide values for the rate constants used to model Ca2+/CaM binding, CaM buffering, Ca2+/CaM/CaMKII as well as Ca2+/CaM/CaN interactions.
Table 5

Rate constants used to model Ca/CaM binding and CaM buffering

Rate Constant

Value

Rate Constant

Value

k 02 C M

4.8387

k 42 B C M

k 42 C M

k 20 C M

10.0*

k 0 B o n C M

3.5 × 10-4‡

k 24 C M

3.4722

k 0 B o f f C M

1.4 × 10-6‡

k 42 C M

500.0*

k 2 B o n C M

k 0 B o n C M

k 02 B C M

k 02 C M

k 2 B o f C M

k 0 B o f f C M

k 20 B C M

k 20 C M /100

k 4 B o n C M

k 0 B o n C M

k 24 B C M

k 24 C M

k 4 B o f f C M

k 0 B o f f C M

Table 5: Rate constants used to model Ca/CaM binding and CaM buffering. Adopted (*) or estimated (‡) from Saucerman et. al. [38].

Table 6

Rate constants used to model Ca/CaM/CaMKII interactions

Rate Constant

Value

Rate Constant

Value

k PP 1

1.72

k 32 C K

2.2

k mPP 1

11.5

k 45 C K

3.35 × 10-3

k 21 C K

65.67164

k 54 C K

3.4722

k 12 C K

328.3582

k 46 C K

2.2 × 10-3

k 13 C K

3.4722

k 64 C K

65.67164

k 31 C K

3.35

k 56 C K

328.3582

k 23 C K

65.67164

k 65 C K

65.67164

Table 6: Rate constants used to model Ca/CaM/CaMKII interactions. All the above parameters are adopted from Saucerman et. al. [38].

Table 7

Rate constants used to model Ca/CaM/CaN interactions

Rate Constant

Value

Rate Constant

Value

k C a o n C N

2.0

k 0 o n C N

46.0

k C a o f f C N

1.0

k 0 o f f C N

537.966

k 02 C N

4.8387

k 2 o n C N

46.0

k 20 C N

0.0606

k 2 o f f C N

3.2604

k 24 C N

3.4722

k 4 o n C N

46.0

k 42 C N

0.199362

k 4 o f f C N

1.3 × 10-3

Table 7: Rate constants used to model Ca/CaM/CaN interactions. All the above parameters are adopted from Saucerman et. al. [38].

Our previous study [30] 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 ( k 25 d h p r ), and the Ca4CaM dependent inactivation rate function ( k 24 d h p r ). The specific formulas for the Markovian state equations and the modified k 25 d h p r and k 24 d h p r 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 k 36 d h p r and k 63 d h p r 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. [39], and Post and Langer [40] 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 [41].

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 [42] with 4 Ca2+ binding for activation and 3 Ca2+ 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 four-state representations during the simulation as Ca2+ 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 [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 [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 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 [4750] is modeled via the rate functions k 12 r y r , k 41 r y r , k 43 r y r and k 32 r y r ; (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 [51]; 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. [52] and Zucchi et al. [53]. Our initial studies indicated that the repriming rate ( k 32 r y r ) in the RyR gating scheme of Stern et al. [46] 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 k 32 r y r tends to reactivate the channel, resulting in saturated Ca2+ release. Therefore, we utilized a value of k 32 r y r 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 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 [5456]. 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 [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 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. [25] 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 [25] 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 [6567], 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 [66]. This release termination mechanism is incorporated in our combined RyR-Luminal sensor model (Figure 4B and Appendix A1, Equations 73-92).

Due to the lack of in-vivo measurements on individual state transitions, the rate constants in our novel luminal sensor model (Table 8) are chosen such that (i) the triadin mediated sensitization of the RyR channel (via A 1 ls ) provides adequate peak RyR release which translates into the upstroke velocity of the cytosolic Ca2+ transient; (ii) the luminal sensor mediated RyR channel inactivation (via A 1 ls ) causes timely release termination resulting in a cytosolic Ca2+ transient duration, physiological for a rat ventricular myocyte; (iii) the rate of post-release RyR recovery results in appropriate channel refractory characteristics. In the case of channels where there is an interaction between ion flow (which is not at equilibrium) and its gating mechanism the microscopic reversibility criteria does not hold true [68, 69]. The RyR channel which is solely modulated by [Ca2+] (both on the dyadic side ([Ca2+] ryr ) as well as the luminal side ([Ca2+] jSR )) experiences a strong interaction of Ca2+ flow through itself and its gating mechanism described as [Ca2+] ryr induced self-inhibition and the luminal sensor induced inactivation. Hence, the rate constants of the luminal sensor model (Figure 4B-i) are not constrained by microscopic reversibility criteria. However, a stability constraint in the form of, the sum of probabilities of all possible states corresponding to triadin (A 1 ls , I 2 ls , I 3 ls , I 4 ls ) and CS (I 2 ls , I 3 ls , B 5 ls , B 6 ls ) being equal to one is explicitly imposed Appendix A1, Equations 87-88) on the model of the luminal sensor in order to avoid run-off.
Table 8

Parameters used in the luminal sensor model

Rate Constant

Value

Rate Constant

Value

k 12 l s

88.16

k 23 l s

57.9

k 21 l s

4.1

k 32 l s

2.42

k 14 l s

0.5

k 35 l s

150.3

k 41 l s

85.7

k 53 l s

25.5

k 42 l s

2.98

k 52 l s

88.16

k 43 l s

25.5

k 56 l s

1.2

k 34 l s

150.3

k 65 l s

401.7

Table 8: Static rate constants used in the model for the luminal sensor

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 [22], whereas those for Ca2+ binding to troponin were taken from Potter and Zott [72]. Rate constants used to describe Ca2+ binding to calsequestrin were based on the study of Cannell and Allen [73], whereas those for Ca2+ 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 Ca2+ 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].

Ca2+-Uptake

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 [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 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 [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 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. [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].

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+ [83]. 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 [84] 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 [85] 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.

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 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 [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 Ca2+ entering through the I Ca,L channel. The characteristics of this channel are hence examined in detail to understand its voltage and Ca2+ dependent behavior. The ability of the I Ca,L channel to facilitate graded RyR release is noted. The high gain associated with RyR release is quantified. This is followed by a study of the properties of the RyR channel with a particular emphasis on its inactivation mechanism. The Ca2+ mediated interaction between these channels is investigated in detail. This is followed by an effort to understand the role of [Ca2+] jSR in CICR. In particular, the relationship of the peak cytosolic calcium transient and the SR Ca2+ content is examined. The effect of modulatory agents like caffeine and thapsigargin on SR release/uptake is studied to understand the significance of these mechanisms in facilitating a normal CICR.

L-type Ca2+ current (ICa,L)

The L-type DHP-sensitive Ca2+ channel has a key role in initiating CICR. Hence a thorough analysis of its activation and inactivation mechanisms is considered here. Though the activation of the I Ca,L channel is solely voltage dependent once activated, the inactivation of the channel is influenced not only by the trans-membrane voltage but also the Ca2+ concentration in the vicinity of the channel [95]. The relative con-tribution of the voltage-dependent and Ca2+-dependent inactivation pathways in the I Ca,L channel model can be studied by selectively blocking each pathway in the model. Figure 7 shows the model-generated waveform for I Ca,L under three conditions; (i) control case with both voltage and calcium dependent in-activation pathways intact; (ii) ryanodine applied to allow only the voltage-dependent inactivation and the Ca2+-dependent inactivation via I Ca,L self-inhibition in the absence of RyR release; and (iii) Ca2+ substituted with Ba2+ to facilitate only the voltage-dependent inactivation pathway and block all Ca2+ dependent inactivation pathways (k 24 dhpr set to zero). The protocol used in case (i) was a 50 ms pulse with clamp voltages ranging from 10 mv to 40 mv in steps of 10 mv from a holding potential of -40 mV. In case (ii) and (iii) the pulse duration was increased to 200 ms to replicate the data.
Figure 7

I Ca,L & [Ca 2+ ] myo - comparison of model generated and experimental data in the positive voltage range. Figure 7: Model generated and experimentally obtained data from a rat ventricular myocyte clamped to voltages between 10 and 40 mv in steps of 10 mv. The holding potential used is -40 mv. Data is obtained from Dr. Philip. T. Palade's lab. (A) ICa,Lin the presence of Ca2+ release from the Ry-sensitive channel. The pulse duration used is 50 ms. (B) ICa,Lwith Ryanodine administered to block RyR release. The pulse duration used is 200 ms. (C) Global cytosolic calcium transients [Ca2+] myo in the presence of RyR release. The pulse duration used is 50 ms (D) ICa,Lusing Barium as the charge carrier in the absence of RyR. The pulse duration used is 200 ms. Solid traces are the model-generated results and dotted lines represent 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 Ca2+-dependent facilitation (CDF). According to Bers et al. [33] 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 [96], suggesting a critical role for Ca2+-dependent facilitation. Our model generated results also indicate that a stronger CDF (in Figure 7A) in the presence of RyR Ca2+ release causes an enhancement in I Ca,L peak amplitude (compare the model fits to data in Figure 7A (presence of RyR release) with Figure 7B (no RyR release) and Figure 7C (Ca2+ substituted by Ba2+)). It is important to note that the peak of the I Ca,L current is reached after the RyR open probability reaches its maximum value owing to the faster dynamics of the RyR channel. Our simulations indicate that following the inward current peak, Ca2+-dependent inactivation (CDI) is much faster and dominates the response for the duration of voltage clamp (30 ms <t <80 ms). In addition, a comparison of Figs. 7A and 7B shows that, Ca2+-dependent inactivation (CDI) caused by Ca2+ release from Ry-sensitive channels is much more significant than the self-inhibition produced by Ca2+ influx via I Ca,L itself. The graded behavior of the cytosolic Ca2+ transient is shown in Figure 7C. Voltage dependent inactivation (VDI) is relatively slow compared with CDI and is best seen in Figure 7D where all Ca2+ inactivation effects have been blocked by Ba2+ substitution for Ca2+. Under RyR blockade (Figure 7B), the relatively slow VDI has its major effect during the late phase of the long voltage clamp pulse (e.g., beyond 100 ms) and in a time range where CDI is relatively constant. Quantitative analysis of movement of states (during the depolarizing pulse duration of 50ms) via different pathways in the six state Markovian model shows that 68.24% of the total inactivation of I Ca,L,TT is via the calcium dependent O 2 dhpr -C 4 dhpr pathway, owing to the large Ca2+ concentration which the channel is exposed to in the dyad. In contrast, only 19.21% of the total inactivation of I Ca,L,SL is via the calcium dependent O 2 dhpr -C 4 dhpr pathway, because of the low cytoplasmic Ca2+ concentration in the vicinity of the channel. Considering the total I Ca,L current (I Ca,L,TT + I Ca,L,SL ), around 63.34% of the I Ca,L channel inactivation is via the Ca2+ dependent pathway.

Figure 8A shows a comparison of model generated and experimentally obtained I Ca,L data (obtained from Dr. Palade's lab and pre-processed to eliminate transients produced by changing clamp voltage in order to obtain model fits using a non-linear least-squares method [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 Ca2+ transient is shown in Figure 8C. Figure 8B shows the well known [52] bell-shaped dependence of the peak [Ca2+] myo on the clamp voltage.
Figure 8

I Ca,L & [Ca 2+ ] myo - comparison of model generated and experimental data in the negative voltage range. Figure 8: Model generated and experimentally obtained data from a rat ventricular myocyte clamped to voltages between -30 and 0 mv in steps of 10 mv. The holding potential used is -40 mv. Data is obtained from Dr. Philip. T. Palade's lab. (A) ICa,Linvoked by negative clamp voltages in the range of -30 mv to 0 mv in steps of 10 mv. The pulse duration used is 100 ms. (B) Peak ICa,Ldependance on clamp voltage. (C) Global cytosolic calcium transients [Ca2+] myo at negative (-30 mv ≤ v ≤ 0 mv) clamp voltages. The pulse duration used is 100 ms (D) Peak cytosolic calcium transients ([Ca2+] myo ) dependence on clamp voltages. Solid traces are the model-generated results and dotted lines represent the data

Figure 9 shows model generated normalized I Ca,L obtained by a voltage clamp to 10 mv from a holding potential of -40 mv using the following protocols. Case 1: Normal release is allowed where, following CICR, a strong Ca2+ dependent inactivation inhibits the I Ca,L current as seen in the plot numbered 1 (data obtained from Dr. Palade's lab). Case 2: Ryanodine applied to block RyR release; Ca2+ entering via the I Ca,L channel causes Ca2+ induced inactivation, although the magnitude of inhibition is far less than in the control case. Case 3: Ba2+ substitution for Ca2+ to completely suppress the Ca2+ inactivation mechanism; the only inactivation pathway present is the slow voltage dependent pathway [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 Ca2+ release (a result of enhanced CDF in the presence of elevated Ca2+ levels). This highlights the important role of Ca2+ in regulating I Ca,L channel opening and thus controlling the amount of extracellular Ca2+ entering the cell.
Figure 9

I Ca,L inactivation pathways. Figure 9: ICa,Linactivation pathways. Voltage clamp protocol used is a 50 ms step pulse to 10 mv from a holding potential of -40 mv. Both model generated results and the data are normalized for easy comparison (1) Control case where both voltage and calcium dependent inactivation pathways are intact. Solid trace is the model-generated result and the dotted line represents the data obtained from Dr. Philip. T. Palade's lab; (2) voltage-dependent inactivation exists however Ca2+-dependent inactivation is only via ICa,Lself-inhibition without RyR release; (3) Only voltage-dependent inactivation pathway preserved by blocking all Ca2+ dependent pathways.

ICa,L,TT - dependent Graded SR-release

One of the most important characteristics of the CICR mechanism is the modulation of the RyR Ca2+ release based on the amount of trigger calcium delivered via the DHPR channel. The voltage-controlled I Ca,L channel is the link between the extracellular excitation and the intracellular Ca2+ release. As shown in Figure 10A, our model reproduces graded release. It is important to note that the onset of the Ca2+ transient is also modulated based on the magnitude of the trigger calcium available to initiate release, as shown in Figure 10A-iii. This is due to the fact that the rate of increase in the open probability of the Ry-sensitive Ca2+ release channel is controlled by the amount of trigger Ca2+ present at the 'mouth' of this channel, which in-turn is graded by the amount of trigger Ca2+ entering the DCU by means of the DHP-sensitive I Ca,L,TT channel. This mechanism is incorporated in our model by allowing the rate constants in the 4-state Markovian model of the Ry-sensitive Ca2+ release channel (Figure 4) to be [Ca2+] ryr dependent. Small changes in i Ca,L,TT cause modulation of pre-release [Ca2+] ryr , which dictates the propensity of the RyR channel for a Ca2+ release. The pre-release values of the rate constants in the 4-state Markovian model of the RyR channel (nonlinear dependence on [Ca2+] ryr as shown in Appendix A1, Equations 79-84) along with the diastolic [Ca2+] jSR (which is kept constant in Figure 10) set the peak RyR open probability achieved by the RyR channel as well as the peak [Ca2+] myo . The use-dependent adaptation of the RyR channel [13] is reflected in its non-linear response to trigger Ca2+. As shown in Figures 10B-iii, iv, the delay in the onset of cytosolic Ca2+ transient closely follows the modulation of the onset of RyR release. Though occurring in separate compartments, the peak of cytosolic Ca2+ transient also tracks the maximum value attained by the open probability of the Ry-sensitive Ca2+ release channel. Here, the clamp voltage was held constant at 10 mv to avoid the interference of voltage dependent Ca2+ transport via the Na+/Ca2+ exchanger which is co-located [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.
Figure 10

Graded CICR. Figure 10: SR calcium release graded by the trigger current ICa,L,TT. (A) - (i) Decreasing trigger current ICa,L,TT(ii) With decreasing trigger current, RyR open probability shows a graded decrease in its peak as well as an increasing delay in its onset (iii) Cytosolic calcium transient follows the trend shown by the RyR release channel (a graded decrease in peak and increasing delay in onset) (iv) RyR open probability vs [Ca2+] ryr for increasing trigger current ICa,L,TT. (B)- (i)Concentration profile in the dyad as a result of RyR release caused by a largest trigger current (see part A(i)) (ii)Concentration profile in the dyad as a result of RyR release caused by a smallest trigger current (see part A(i)) (iii) With decreasing trigger current, RyR open probability shows a graded decrease in its peak as well as an increasing delay in its onset (shown on an expanded time scale) (iv)Difference in concentration between the mouth of the DHP and Ry-sensitive Ca2+ channels. (Inset shows the cross-over when this difference becomes negative indicating trigger current dependent latency involved in 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. [13] 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 [13]. 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 [99], a Ca2+ spark translates into a release from around 4 RyR channels (also reported by Blatter et. al. [28]). 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 [25] that the luminal sensor plays a fundamental role in an active extinguishing mechanism [51] that effects a robust [Ca2+] jSR - dependent closure of the RyR channel. This mechanism for inactivation of the Ry-sensitive Ca2+ release channel is accounted for in our model.

Cytosolic Ca2+ transient is also graded by the duration of the I Ca,L,TT trigger current controlled by the voltage clamp pulse duration. Our model also shows that triggered release can be prematurely stopped by rapid repolarization [100102]. The effect of depolarization duration on the time course of the cytosolic Ca2+ transient is indicated in Figure 11 where the duration of the voltage clamp pulse is decreased from 80 ms to 5 ms (while keeping the clamp voltage constant at 10 mv) resulting in premature stoppage of release and thus a decrease in peak cytosolic Ca2+ transient. This effect is a combined result of (i) the modulation of release due to pulse duration dependent change in the amount of trigger Ca2+ entering the dyadic coupling unit (DCU) and (ii) the pulse duration dependent change in the relative role of the Na+/Ca2+ exchanger co-located in the DCU.
Figure 11

Effect of depolarization duration on [ Ca 2+ ] myo . Figure 11: The effect of changing the duration of depolarizing pulses T pd on the time course of the [Ca2+] myo transient. Voltage clamp protocol used is a step pulse to 10 mv from a holding potential of -40 mv. Normalized [Ca2+] myo transient corresponding to depolarizing pulses of 5, 10, 20, and 80 ms duration. This model-generated result shows similarity to data in Figure 2, Cannell et al. [100] showing pulse-duration dependent changes in cytosolic Ca2+ concentration.

High Gain of Ca2+ Release

Besides the graded release, an extremely valuable characteristic of the CICR process is the high gain associ-ated with it. A small amount of Ca2+ entering the DCU via the I Ca,L,TT channel causes a large release of Ca2+ from the sarcoplasmic stores via the Ry-sensitive Ca2+ release channel on the jSR lumen that interacts with the DCU. This high gain is essential in producing physiological cytosolic Ca2+ transients when 10000 of these DCU's operate in tandem. Two different definitions of the gain or amplification factor due to CICR have been adopted in our model simulations, namely the ratio of:
  1. 1.

    average integrated RyR flux to average integrated DHPR flux [7]; and

     
  2. 2.

    the peak Ca2+ 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 Ca2+ flux through DHPR and RyRs, respectively. Figure 12B shows cytosolic calcium transients under conditions of ryanodine blockade and with RyR Ca2+ release. Based on our measurements, there is approximately a 2 ms delay from the onset of the DHPR influx to the initiation of the RyR release flux.
Figure 12

CICR Ca 2+ gain. Figure 12: High gain associated with CICR. Voltage clamp protocol used is a 50 ms step pulse to 10 mv from a holding potential of -40 mv. (A) Gain formulation by comparison of Ca2+ flux through the DHP sensitive and RyR sensitive Ca2+ channels. The voltage dependence of gain measured using this method is shown in the inset. (B) Gain formulation by comparison of cytosolic calcium transient formed in the presence and absence of RyR Ca2+ release.

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 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. [104]). It is important to note that, with increasing clamp voltage, the decreased ability of the Na+/Ca2+ exchanger (which is co-located [98] in the dyad) to extrude Ca2+ partially compensates for the declining trigger current, in facilitating SR release and hence assists in increasing CICR gain.

RyR Refractory Characteristics

Given the fact that the CICR mechanism has a high gain associated with it, it could be prone to unstable behavior if not for the refractory characteristics of the Ry-sensitive release channel [105]. It is now known [25] that the luminal sensor plays a fundamental role in an active extinguishing mechanism [51] effecting a [Ca2+] jSR dependent robust closure of the RyR channel which is accounted for in our model. The RyR release fails to become regenerative due to the role of the luminal sensor which after a release occurs, forces the RyR channel into an absolute refractory state followed by a relative refractory state as shown in Figure 13C. When the RyR channel is in the absolute refractory period, [Ca2+] ryr drops to a level much below what is caused by a sparklet (trigger Ca2+) hence it robustly avoids re-excitation. This extremely critical refractory feature of the channel enabled by the RyR luminal sensor acts as a protective mechanism against premature Ca2+ release in the wake of a secondary excitation that occurs before jSR can be filled back to the control level. The refractory nature of the channel is caused by the [Ca2+] jSR dependent inhibition induced by the protein triadin from the luminal side of the RyR channel. This restraining lock on the RyR channel assists in reloading the jSR.
Figure 13

RyR refractory characteristics. Figure 13: RyR refractory characteristics. (A) The stimulation protocol comprises of two consecutive stimuli S1 & S2, both of which are identical to the dyadic component (IC a,L,TT) of trace 1 (ICa,Levoked by a 10 mv clamp pulse) in Figure 7A with a peak amplitude of 6.657 pA/pF and a duration (T1) of 50 ms. These are applied with their separation in time (T2) increasing from 50ms to 250 ms in steps of 25 ms. (B) The RyR sensitive Ca2+ release channel shows refractory characteristics with an absolute refractory period of 50 ms and a relative refractory period of 250 ms. (C) The corresponding [Ca2+] jSR traces reflect partial RyR recovery. (D) Cytosolic calcium transient recovers completely after 250 ms. The dotted line indicates lack of RyR release when T2 = 75 ms due to insufficient RyR recovery.

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 [106] on rat ventricular myocytes. When the RyR channel is in the absolute refractory period, [Ca2+] RyR drops to a level much below what is caused by a sparklet (trigger Ca2+) and hence robustly avoids re-excitation. However, when the RyR channel is in the relative refractory period, a partial release (Figure 13C) is possible, despite the RyR receptors affinity for release being low.

Figures 14A and 14B show the steady state [Ca2+] myo and [Ca2+] jSR predicted by the model, with and without a functioning luminal sensor, respectively. The stimulation protocol used is a pulse train of amplitude (-40 mv to 10 mv), duration (50 ms) and frequency of 4.0 Hz which is applied for a period of 2.5 sec. The value of the luminal control of the RyR channel in allowing adequate SR filling can be seen in Figure 14 where, in the presence of the luminal sensor the cell maintains a normal cytosolic Ca2+ transient (Figure 14A-i). This is made possible due to sufficient SR filling (Figure 14A-ii) as a result of RyR refractoriness. However, inactivation of the RyR channel is also known to depend on the high local Ca2+ concentration consequential to its own Ca2+ release [51]. The lack of a luminal sensor forces the RyR channel to solely rely on Ca2+ dependent inactivation mechanism. The resulting inadequate RyR inactivation depletes diastolic [Ca2+] jSR level (Figure 14B-ii) causing the cytosolic Ca2+ transient to diminish to new lowered steady state values (shown in Figure 14B-i). The absence of the luminal sensor is modeled by setting the value of the luminal sensing variable ('var' in the model) to a level consistent with it's normal value under 4 Hz, voltage clamp stimulation. This mimics the case where the luminal sensor is insensitive to changes in [Ca2+] jSR and hence non-functional. It is important to note that the decrease (37.07%) in diastolic SR level (as indicated in Figure 14B-ii) is due to the lack of a luminal control on the RyR channel resulting in a reduced rate of RyR recovery which in turn reflects in a compromised SR filling rate. Hence, the presence of a luminal sensor which monitors the SR Ca2+ content, is key to long term Ca2+ stability in the cell.
Figure 14

Role of Luminal sensor. Figure 14: Application of a train of voltage clamp pulses (amplitude -40 to 10 mV; duration 50 ms) for a period of 2.5 sec. Pulse repetition rate is 4.0 Hz. Luminal sensor is intact in A and removed in B. Both A and B are steady state responses to voltage clamp stimulation. A-(i) [Ca2+] myo transients with normal amplitude. A-(ii) Sustained SR filling due to robust RyR inactivation. B-(i) Diminishing [Ca2+] myo transients indicate poor SR filling. B-(ii) Lowered [Ca2+] jSR levels indicating inadequate refilling of SR stores due to lack of RyR refractoriness.

CICR modulation by the jSR Ca2+ content

From the refractory characteristics of the RyR channel, it is evident that the RyR release is strongly dependent on the jSR Ca2+ content. In fact, the SR Ca2+ release through the RyR channel depends on: (a) the amount of trigger Ca2+ entering by means of the I Ca,L,TT channel (b) the concentration gradient of Ca2+ between the jSR and the mouth of the RyR channel in the dyadic space and (c) the open probability of the RyR channel modulated by the interaction between the luminal sensor and the RyR protein. The SR Ca2+ released also inactivates the trigger current I Ca,L,TT (as shown in Figure 7), hence causing an indirect self-inhibition of RyR release.

Figure 15 shows a phase plot of RyR open probability (O 2 ryr ) versus the Ca2+ concentration at the mouth of the RyR channel in the dyadic space [Ca2+] ryr constructed from model-generated data for a pulse of amplitude -40 mv to 10 mv and a duration of 50 ms. The pre-release diastolic [Ca2+] jSR is 918 μ M. Following excitation by the trigger current during phase A, the RyR channel begins to open, allowing Ca2+ flux into the dyad and elevating [Ca2+] ryr at the mouth of the RyR channel. The rate of opening of the channel is Ca2+ dependent; hence as [Ca2+] ryr increases, RyR open probability increases first slowly and then ever more rapidly. With this onrush of Ca2+, [Ca2+] ryr rapidly equilibrates with the [Ca2+] jSR , reducing and eventually abolishing the concentration gradient across the channel. [Ca2+] ryr soon reaches its maximum value (T1 in Figure 15) and begins to track the decrease in [Ca2+] jSR , there being no further significant concentration gradient existing between [Ca2+] jSR and [Ca2+] ryr at the mouth of the RyR channel. During phase B, the rate of rise in RyR open probability decreases due to decreasing [Ca2+] ryr and the onset of self-inhibition due to the large values of Ca2+ concentration at the mouth forcing the RyR open probability to attain its maximum value (T2 in Figure 15). This is followed by phase C, where the RyR open probability begins to decrease slowly, initiating channel recovery due to decreasing overall [Ca2+] dyad levels as a result of (a) Ca2+ fluxing out of the dyad into the cytosol, (b) lack of a drive from SR release due to drastically diminished Ca2+ gradient at the mouth of the RyR channel and (c) continued Ca2+ induced self-inhibition. Decreasing levels of free Ca2+ in the jSR force a strong RyR inactivation via the luminal sensor (beginning at T4 and continuing throughout phase D; Figure 15), which in turn causes a sharp decline in RyR open probability. The minimum value of [Ca2+] jSR occurs at T5, and the channel is ultimately closed at T6 (Figure 15). The luminal sensor mediated inactivation assists in robust [Ca2+] jSR recovery.
Figure 15

Ryanodine receptor open probability. Figure 15: RyR open probability modulation by concentration at the mouth. The pre-release diastolic [Ca2+] 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, [Ca2+] ryr initially increases rapidly with a small accompanying rise in RyR open probability followed by steeper increase to peak [Ca2+] 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 [Ca2+] ryr level and onset of Ca2+ induced self-inhibition; Phase C - Decreasing Ca2+ dyad and continued Ca2+ induced self-inhibition, begins the channel recovery process; Phase D - Decreasing [Ca2+] jSR forces a strong RyR inactivation via the luminal sensor hence robustly closing the channel.

Cytosolic peak [Ca2+] dependence on SR Ca2+ content

Figure 16A shows three characteristic regions in the plot of the peak [Ca2+] myo vs [Ca2+] jSR . In Region I, where the jSR load is small, two things occur: (a) RyR release is reduced, as is the basal intracellular Ca2+ concentration, and (b) Ca2+ dependent inactivation of the I Ca,L,TT channel is reduced due to the lowered dyadic Ca2+ concentration, which allows for greater Ca2+ entry into the dyadic space via the I Ca,L channel. This increase in trigger current results only in a very small increase in Ca2+ release, because of the inhibition of the Ry-sensitive channel caused by the luminal sensor in response to low SR Ca2+ content. Thus, these factors cumulatively result in only a small linear increase in the peak [Ca2+] myo with increasing SR Ca2+ content in Region I. The behavior in Region I corresponds to traces 7-1 in Figure 16B (shows phase plots of O 2 ryr vs [Ca2+] ryr for different steady state [Ca2+] jSR concentrations), where the peak RyR open probability shows very small gradual increase with increasing levels of diastolic pre-release [Ca2+] jSR for a constant trigger Ca2+ input via I Ca,L . The steady state diastolic [Ca2+] jSR levels being low, the maximum RyR open probability values attained are low. Although the peak [Ca2+] ryr does not increase substantially as the diastolic [Ca2+] jSR levels increase from 317μ M to 673μ M, the area enclosed by the [Ca2+] ryr - O 2 ryr loop (which indicates the amount of SR Ca2+ released into the dyad) increases exponentially due to increasing peak RyR open probability combined with the increasing difference between rate of activation and inactivation with rate of activation increasing faster than the rate of channel inactivation. In Region II of Figure 16A, increasing SR Ca2+ content begins to translate into a substantial increase in peak [Ca2+] myo . This is because increasing diastolic [Ca2+] jSR causes a significant increase in the peak RyR open probability, as shown in Figure 16C, which translates into a commensurate increase in the peak of the cytosolic Ca2+ transient. The inhibiting role of the luminal sensor is relieved with building [Ca2+] jSR levels. The nonlinearity observed between the diastolic [Ca2+] jSR levels 673μ M (Figure 16B) and 734μ M (Figure 16C) at the transition between Regions I and II corresponds to the existence of a threshold characteristic for RyR release [107, 108]. As the diastolic [Ca2+] jSR levels increase beyond 734μ M, the area enclosed by the [Ca2+] ryr - O 2 ryr loop increases substantially reaching a maximum value at [Ca2+] jSR level of 948μ M not only due to a significant increase in both peak [Ca2+] ryr and peak O 2 ryr but also due to a rapid increase in rate of RyR activation resulting in a much larger difference in rate of activation and inactivation. Region III exhibits a large rapid increase in peak [Ca2+] myo due to the large RyR Ca2+ release associated with the high SR Ca2+ content. However, it is important to note that this is not a result of increasing peak RyR open probability, which begins to plateau for steady state [Ca2+] jSR values larger than 850μ M (traces 6-1 in Figure 16C). The continuing rapid increase in peak cytosolic Ca2+ concentration in region III is a combined result of: (a) increased release owing to large SR Ca2+ content; (b) large values for saturated RyR open probability supported by an increase in the area contained by each loop as seen in traces 6-1 of Figure 16C; and (c) saturated operation of the SERCA pump which acts as a predominant buffer in restoring the Ca2+ concentration levels in the cytosol after RyR release [77]. As the diastolic [Ca2+] jSR levels increase beyond 948μ M, the area enclosed by the [Ca2+] ryr - O 2 ryr loop begins to decrease despite increasing peak [Ca2+] ryr due to saturation of the peak RyR open probability and increasing rate of RyR channel inactivation due to large values of local Ca2+ concentration ([Ca2+] ryr ) assisting in faster recovery [51]. This model generated relationship between the peak [Ca2+] myo and the pre-release [Ca2+] jSR agrees with Figure 4 in Diaz et. al. [29].
Figure 16

Peak [Ca 2+ ] myo vs [Ca 2+ ] jSR . Figure 16: (A) Relation between the peak of the cytosolic Ca2+ transient and the diastolic [Ca2+] jSR . Region I: The peak RyR open probability is very small owing to the low SR Ca2+ content. This causes the peak Ca2+ transient to be very small. Region II: Increasing SR Ca2+ content causes an increase in the peak RyR open probability and hence in the peak of the cytosolic Ca2+ transient. Region III: A large rapid increase in the peak Ca2+ myo is observed due to the large RyR Ca2+ release associated with the high SR Ca2+ content and saturation of RyR open probability. This model-generated result shows similarity to data in Figure 4, Diaz et al. [29].(B) Dependence of RyR open probability (O2 RyR ) on diastolic Ca2+ jSR . Increasing pre-release diastolic Ca2+ jSR results in increasing peak O2 RyR . However, owing to the low SR Ca2+ levels the peak open probability is very small. The stimulation protocol used is a pulse of amplitude -40 mv to 10 mv with a duration of 50 ms. (C) Dependence of RyR open probability (O2 RyR ) on diastolic Ca2+ jSR . Increasing diastolic Ca2+ jSR results in a linearly increasing peak O2 RyR . In the open probability saturation region the increase in peak cytosolic Ca2+ transient is sustained by increasing the area contained by each loop as seen in traces numbered 5 to 1. The progress in time occurs along the arrow indicated on trace 1. The stimulation protocol used is a pulse of amplitude -40 mv to 10 mv with a duration of 50 ms.

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

Caffeine

Adachi-Akahane et al. [71] investigated I Ca,L -induced Ca2+ release in rat ventricular myocytes in the presence and absence of caffeine, which altered [Ca2+] myo . Using their protocols and data as a guide, we modeled the effect of caffeine on jSR Ca2+ release and its subsequent effect on I Ca,L . Specifically, a Boltzman relationship was adopted to simulate the binding of caffeine to the ryanodine receptor increasing channel conductance, which in turn enhances jSR Ca2+ release. The process is assumed to be dose-dependent, with an open probability of 0.5 at a caffeine concentration of 5 mM (Appendix A1, Eq. 86).

Figures 17A and 17B show model-generated [Ca2+] myo responses in myocytes dialyzed with 2 mM Fura-2, where a simple voltage clamp pulse (-40 mV to 0 mV) is applied in the presence and absence of caffeine. I Ca,L waveforms associated with these two protocols are shown in Figure 17C. In protocol A, conducted at normal rest levels of [Ca2+] myo (10 nM), the depolarizing test pulse to 0 mV from the holding potential of -40 mV for 0.1 sec fully activates I Ca,L , which in turn triggers a rapid but small amplitude Ca2+ transient (a rise from 30 nM to 95 nM as shown in Figure 17A. Protocol B starts with an application of 5 mM caffeine for a period of 1.0 sec (open state probability of RyR channel is 0.5), followed by the same depolarizing test pulse applied at 1.0 sec. The conditioning caffeine pulse induces a [Ca2+] myo transient rising from a resting value of 30 nM to a peak value of 110 nM (Figure 17B, first peak -note the different timescale). The short depolarizing pulse activates I Ca,L , which triggers a much smaller [Ca2+] myo transient via CICR (30 nM, Figure 17B, second peak). This I Ca,L -evoked Ca2+ transient in response to the voltage pulse is smaller in amplitude than control (Figure 17A), due to: (a) reduction in RyR release due to depletion in [Ca2+] jSR ([Ca2+] jSR drops from 1.25 mM to 0.25 mM with caffeine application).; and (b) the strong competing role of the luminal RyR sensor in inactivating the RyR channel in response to the post-caffeine, low SR content. The I Ca,L -inactivation characteristics undergo a change as shown in Figure 17C, where trace B inactivates more slowly compared to the control I Ca,L (trace 1). After exposure to caffeine, the overall amplitude of [Ca2+] dhp during release does not increase as much as in the control case, since jSR Ca2+ release is reduced. Consequently, Ca2+ dependent inactivation (during release) of I Ca,L after caffeine application is not as strong as that in the control case. However, the baseline level of Ca2+ concentration at the mouth of the DHP channel ([Ca2+] dhp ) is elevated following caffeine application, causing sustained inactivation resulting in the crossover of the current trace in Figure 17C. This is consistent with the observation of Cens et al. [41], that Ca2+ induced inactivation is a result of calmodulin mediated sensing of the local Ca2+ concentration.
Figure 17

Effects of the application of caffeine. Figure 17: Model-simulated effects of the application of caffeine on ICa,L,TTand [Ca2+] i in a cell dialyzed with 2 mM Fura-2. Changes in [Ca2+] 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].

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 Ca2+ uptake is blocked, SR Ca2+ content declines, as does the jSR Ca2+ release with each voltage pulse (Figure 18A). Considerable time is consumed after thapsigargin application, before SR Ca2+ content is finally exhausted (8-10 beats at 0.05 Hz; Figure 18A). As the magnitude of the [Ca2+] myo transient decreases, the [Ca2+] dhp dependent inactivation of I Ca,L slows (Figure 18B), resulting in an enhanced peak.
Figure 18

Effects of the application of Thapsigargin. Figure 18: In rat ventricular myocytes dialyzed with 2 mM Fura-2, thapsigargin gradually abolishes [Ca2+] myo transients and decreases the rate of inactivation of I Ca,L . This result is generated by the model, in response to a pulse train of amplitude (-40 mV to 0 mV), duration (0.1 sec) and frequency of 0.05 Hz. (A) Time course of the effect of thapsigargin on Δ[Ca2+]i. (B) L-type Ca2+ currents at the times indicated in (A). This model-generated result shows similarity to data in Figure 11, Adachi-Akahane et al. [71]

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 [109].

Effect of modulation of [Ca2+]o

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

Effects of modulating [Ca 2+ ] o . Figure 19: The effects of changing external Ca2+ concentration on cytoplasmic Ca2+ and SR Ca2+ 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 [Ca2+] o the instantaneous increase in [Ca2+] myo occurs due to more Ca entry, followed by a gradual staircase rise due to the the change in SR content; (B) Changes in [Ca2+] myo forced by a sudden decrease in [Ca2+] o ; (C) Changes in [Ca2+] jSR due to increased [Ca2+] o ; (D) Changes in [Ca2+] jSR due to reduced [Ca2+] o . This model-generated result shows similarity to data in Figure 2, Diaz et al. [110].

Calcium balance under conditions of repetitive stimulation

Our model exhibits long term Ca2+ 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 steady-state operation, Ca2+ entry into the cytosol via I Ca,L must exactly balance Ca2+ efflux. Changing the rate or pattern of stimulation can have significant effects on the cell's cytosolic Ca2+ balance and subsequent contractile response [16]. To demonstrate Ca2+ balance in our model, long (Figure 20A) or short (Figure 20B, C and 20D) voltage clamp pulses were applied at a selected repetition frequency, and the resultant changes in [Ca2+] myo were correlated with sarcolemmal and SR Ca2+ fluxes. Analogous to the protocol of Negretti et al. [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 (short-pulse) at a frequency of 0.33 Hz, for a total test duration of 2 minutes (40 pulses). In analyzing the Ca2+ balance, the amount of electric charge Q transferred per pulse was calculated according to the equation: 0 T I ( t ) d t , where I(t) can be I Ca,L , I NaCa , I PMCA , I ryr , I cyt,serca , and T is 1 cycle duration (3 ms for 0.33 Hz stimulation).
Figure 20

Long term Ca 2+ stability during repetitive stimulation. Figure 20: The stimulation protocol used is a pulse train of amplitude (-40 mv to 10 mv) and frequency of 0.33 Hz.(A) Integrals of external and internal Ca2+ currents on stimulation with repetitive long pulses (800 ms); (B) Percentage of changes in [Ca2+] myo peak magnitude and cumulative SR Ca2+ loss on stimulation with repetitive short pulses (100 ms); (C) Integrals of external Ca2+ currents with short-pulse stimulations (100 ms); (D) Integrals of internal Ca2+ currents with short-pulse stimulations (100 ms). This model-generated result shows similarity to data in Figure 3 and Figure 5, Negretti et al. [111].

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+).

Long-pulse protocol

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.

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 [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. [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 Ca2+ concentration [Ca2+] jSR . Both of these indicators show a net loss of Ca2+ from the cell. Figure 20C and 20D show the transient charge movements associated with the external and internal Ca2+ currents respectively during application of the short-pulse train, which are indicative of an elevated Ca2+ load and resultant Ca2+ imbalance when the pulse train is first applied. Balance is achieved, but at new lower values of [Ca2+] myo and [Ca2+] jSR . For convenient comparison, all internal fluxes are made positive and all external fluxes negative. Figure 20C shows that I NaCa plays a significant role in forming the transient response seen in the short-pulse protocol in extrusion of excess Ca2+ load on the cell. The decreasing amplitude of the Ca2+ transient causes a decreased LSR uptake current (I cyt,serca ), which in turn causes a decrease in [Ca2+] jSR resulting in decreased release (I ryr ). With successive pulses, I cyt,serca and I ryr decline further, mirroring the exponential decline in I NaCa toward balance. At steady state, the external and internal fluxes reach the new equilibrium, where the sum of the integrals of external currents and the sum of the integrals of internal currents are zero.

Long-term calcium stability at higher pacing rates

Most VC experiments on cardiac cells are usually conducted at low stimulation rates. Our simulations of typical VC experiments at these low rates all exhibit long term calcium stability (Figures 17, 18, 19 and 20). However, our model can also exhibit sustained calcium balance at higher pacing rates. Figure 14A shows steady state cytosolic Ca2+ transients in response to a repetitive 4 Hz (which is more physiological for a rat ventricular myocyte) voltage clamp stimulation lasting for 2.5 seconds. The sustained peak and basal Ca2+ levels of [Ca2+] myo over this prolonged time frame are indicative of long tern Ca2+ stability in the model. The corresponding [Ca2+] jSR profile over 2.5 seconds is also shown in Figure 14A, which indicates sustained SR filling owing to the luminal sensor mediated control of the RyR channel. The stimulation protocol used is a pulse train of amplitude (-40 mv to 10 mv), duration (50 ms) and frequency of 4.0 Hz. Tables 9 and 10 provide values for the initial conditions used in our model.
Table 9

Initial Conditions

State variable

Definition

Initial Value

[Ca2+] myo

Ca2+ concentration in the myoplasm

8.1027 × 10-5 mM

[Ca2+] jSR

Ca2+ concentration in jSR

1.2677 mM

[Ca2+] LSR

Ca2+ concentration in LSR

1.3346 mM

[Na+] myo

Na+ concentration in myoplasm

16.746 mM

[Na+] dyad

Na+ concentration in dyadic space

16.321 mM

[Cs+] myo

Cs+ concentration in myoplasm

140.2154 mM

[Cs+] dyad

Cs+ concentration in dyadic space

140.2157 mM

O c

Fractional occupancy of Calmodulin by Ca2+

0.033091

O tc

Fractional occupancy of Troponin by Ca2+

0.016049

O tmgc

Fractional occupancy of Troponin by Mg2+ and Ca2+

0.321764

O tmgmg

Fractional occupancy of Troponin by Mg2+

0.598385

[CaF 3]

Concentration of Fluo3-Ca2+ complex

21.88721 μM

C 1 ryr

Closed (resting) state of RyR channel

0.9990953

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

0.1673614

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

0.6709816

I 2 ls

Fraction of Tr/J bound to RyR and Calsequestrin

0.155258

I 3 ls

Fraction of Tr/J bound only to Calsequestrin

0.0623695

B 6 ls

Fraction Calsequestrin bound to Ca2+

0.619346

[Ca2+] serca

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

PLB dp

Fraction of unphosphorylated Phospholamban

7.684160 × 10-2

Ca 2 CaM

2 Ca2+ ions bound to C-terminus of CaM

34.56529

Ca 4 CaM

4 Ca2+ ions bound to C & N terminus of CaM

8.635052 × 10-2

CaMB

Buffered CaM

7.563836 × 10-2

Ca 2 CaMB

Buffered Ca2CaM

2.035086

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

CaMCaN

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

3.473412μM

[Ca2+] dyad

Ca2+ concentration in the dyadic space

9.012 × 10-5 mM

 

Table 9: Initial conditions used for the state vector in order to solve the linearized system of ordinary differential equations

Table 10

Initial Conditions

State variable

Definition

Initial Value

P 1

Fraction of inactive dephosphorylated CaMKII in Ca2CaM bound state

5.527608 × 10-1

P 3

Fraction of active dephosphorylated CaMKII in Ca4CaM bound state

3.661260 × 10-1

P 6

Fraction of active Thr287-autophosphorylated but CaM autonomous CaMKII

1.314410 × 10-3

P 5

Fraction of active Thr287-autophosphorylated but Ca2CaM bound CaMKII

6.277911 × 10-7

P 4

Fraction of active Thr287-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

Secondary [Ca2+]myo transients induced by "tail currents"

Our voltage clamp simulations have dealt with clamp voltages in the range -30 mV ≤ V ≤ 40 mV. Some experimental studies [100, 2, 112] however have employed even larger clamp voltages (40 mV ≤ V ≤ 60 mV) to explore the high voltage behavior of the DHP-sensitive calcium channel. Such large clamp voltages elicit a brief Ca2+ influx called a "tail current", which has been shown to trigger RyR release and hence cause contraction during repolarization [100, 2]. The secondary Ca2+ transient induced by the "tail current" is a critical argument in favor of CICR as these "tail currents" are not observed in skeletal muscle where the membrane potential directly controls SR Ca2+ release [112]. Our model reproduces this secondary Ca2+ transient observed during the return to resting potential from a large clamp voltage (≥ 40 mV). Figure 21A shows a cartoon depicting the voltage clamp stimulation protocol used in our study where a pulse of amplitude (-40 mV to +50 mV) is employed with the pulse duration (T p ; Figure 21A) increasing from 50 ms (trace 1) to 200 ms (trace 7) in steps of 25 ms. The clamp potential transition time (T t ) is fixed at 1 ms. The peak of the "tail current" at the end of the pulse is 25 times larger than the peak of the I Ca,L current observed at the beginning of the pulse. This is in agreement with model generated VC data reported by Geenstein et al.([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 Ca2+ channel, with RyR Ca2+ release indicated by the corresponding decrease in [Ca2+] jSR (panel C) and increase in [Ca2+] myo (panel D). We note from Figure 21C that for a shorter 50 ms VC pulse, there is no RyR Ca2+ release although an I Ca,L tail current is produced. This tail current produced with the 50 ms pulse slightly augments the [Ca2+] myo transient just after its peak (trace 1, panel E), but there is no corresponding decrease in the [Ca2+] jSR transient (trace 1, panel D) indicating no RyR Ca2+ release. At the end of this short pulse (50 ms), the RyR channel is in fact in its absolute refractory period. The small secondary increase in [Ca2+] myo transient as seen in data corresponding to trace 4 of Figure 7C is a result of the large Ca2+ influx via the I Ca,L "tail current". During this phase, the Na+/Ca2+ exchanger is biased to extrude Ca2+ out of the dyad, and hence cannot be responsible for this secondary rise in the cytosolic Ca2+ transient. As pulse duration is increased beyond 100 ms the "tail current" causes a gradual increase in RyR open probability (traces 3-7 in Figure 21) indicating progressive recovery from the refractory period. Minimum value of [Ca2+] jSR in traces 3-7 of Figure 21D show increase in SR Ca2+ release with increasing pulse duration as a result of adequate RyR channel recovery and increasing I Ca,L "tail current". Cytosolic Ca2+ transients shown in Figure 21E indicate "tail current" induced secondary RyR release for pulse duration ≥ 100 ms. It has also been previously reported [113] that recruitment of an additional population of previously 'silent' Ca2+ channels could cause facilitation of tail currents at increasingly large clamp voltages with a time-dependence associated with the recruitment process. Our model not only accounts for the contribution of the already open channels (in states O 2 dhpr & O 3 dhpr ) to the tail current during repolarization but also allows for the voltage- and time-dependent recruitment of 'silent' Ca2+ channels modeled using a high voltage state C 6 dhpr . This study proves to be an adjunct to the study presented in Figure 13 in establishing the refractory nature of the RyR channel.
Figure 21

Secondary [Ca 2+ ] myo transients induced by "tail currents". Figure 21: The stimulation protocol used is a pulse of amplitude (-40 mv to +50 mv) with pulse duration T p increasing from 50 ms (trace 1) to 200 ms (trace 7) in steps of 25 ms. (A) A cartoon depicting the voltage clamp stimulus protocol where, pulse duration (T p ) is varied to obtain traces 1-7 and the clamp potential transition time (T t ) fixed at 1 ms. (B) I Ca,L tail currents elicited during repolarization from +50 mv to -40 mv. (C) Open probability of the RyR channel indicating gradual recovery from the refractory period (D) [Ca2+] jSR traces show increase in RyR release with increasing pulse duration as a result of adequate RyR channel recovery. (E) Cytosolic calcium transients indicate secondary RyR release caused by tail currents.

Cytosolic Buffering

In our model, cytosolic buffering is attributed to several factors including: (a) calmodulin (CaM); (b) the Ca2+-specific (Tc) troponin binding site; (c) the Ca2+ - Mg2+ competitive troponin binding site; and (d) the fluorescent indicator dye Fluo3 used to detect changes in [Ca2+] myo . The effects of the component buffers in helping to maintain a low [Ca2+] myo are shown in Figure 22 in terms of occupancy functions, such as O C (fractional occupancy of calmodulin by Ca2+), O tc (fractional occupancy of troponin-Ca sites by Ca2+), O tmgc (fractional ocupancy of troponin-Mg sites by Ca2+), O tmgmg (fractional occupancy of troponin-Mg sites by Mg2+) in the cytosol, and O calse (fractional occupancy of calsequestrin by Ca2+) in the jSR Ca2+ release compartment. Figure 22 shows that when the Ry-sensitive Ca2+ release channel is triggered, the jSR releases Ca2+ and the occupancy O calse in the jSR, declines from 68% to 37%. Ca2+ release from the jSR induces fast Ca2+ binding by calmodulin and troponin in the cytosol, as represented by the increases of O c (0.29) and O tc (0.16) in Figure 22. Interactions between Ca2+ and the troponin-Mg sites result in an increase in occupancy of these sites by Ca2+ (O tmgc ), and a decrease in occupancy by Mg2+ (O tmgmg ).
Figure 22

Occupancy of Ca 2+ buffers. Figure 22: Contributions of fractional occupancy of Ca2+ buffers in SR and cytoplasm. Transfer from O c and O tc to O tmgc occurs as a result of the pulse stimulation. O tmgmg is shown as a dotted line as it does not reflect Ca2+ binding. The corresponding [Ca2+] myo is overlayed (note the different axis)

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 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 [53]. 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. [71] and Diaz et al. [110]). 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 [114].

Model Limitations

  1. (a)

    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.

     
  2. (b)

    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).

     
  3. (c)

    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. ([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 Ca2+ 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 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.

Appendix

Below is the complete set of equations used in the model.

A1 - Equations for currents in the model

L-Type Ca 2+ current

Ca2+ current through the DHP-sensitive ICa,Lchannel
I C a , L = R C a , L ( O 2 d h p r + O 3 d h p r ) P C a Z C a 2 F 2 V R T × ( [ C a 2 + ] d h p e 2 F V R T 341.0 [ C a 2 + ] 0 e 2 F V R T 1 )
(1)
Na+ current through the DHP-sensitive ICa,Lchannel
I N a = 1.056 × e ( ν 21.73 21.23841 ) 1.056 × e ( ν 21.73 21.23841 ) + [ C a o 2 + ] P N a Z N a 2 F 2 V R T × ( 0.75 [ N a + ] d h p e F V R T 0.75 [ N a + ] 0 e F V R T 1 )
(2)
Cs+ current through the DHP-sensitive ICa,Lchannel
I C s = 1.056 × e ( ν 21.73 21.23841 ) 1.056 × e ( ν 21.73 21.23841 ) + [ C a o 2 + ] P C s Z C s 2 F 2 V R T × ( 0.75 [ C s + ] d h p e F V R T 0.75 [ C s + ] 0 e F V k C s R T 1 )
(3)

where k Cs = 0.5 and [Ca2+] dhp , [Na+] dhp , [Cs+] dhp are Concentrations at the mouth of the DHP-sensitive Ca2+ channel.

The corresponding unitary currents i Ca , L , i Na , i Cs are obtained by dividing the above net channel currents by the number of dyadic units, N dyad .

Gating scheme for the 2-state Markovian model used to allow Ca2+ mediated interaction of the DHP-sensitive I Ca , L channel and calmodilin:
S 1 d h p r = 1 S 2 d h p r
(4)
d S 2 d h p r d t = k s 12 d h p r S 1 d h p r k s 21 d h p r S 2 d h p r
(5)
Expressions for the rate constants:
k s 12 d h p r = 4.0 × e ( [ C a 4 C a M ] d h p 100.0 90.0 )
(6)
k s 21 d h p r = 18.0
(7)
Gating scheme for the 6-state Markovian model for the DHP-sensitive L-type Ca2+ release channel:
C 5 d h p r = 1 ( C 1 d h p r + O 2 d h p r + O 3 d h p r ) ( C 4 d h p r )
(8)
d C 1 d h p r d t = k 21 d h p r O 2 d h p r + k 51 d h p r C 5 d h p r ( k 12 d h p r + k 15 d h p r ) C 5 d h p r
(9)
d O 2 d h p r d t = k 12 d h p r C 1 d h p r + k 42 d h p r C 4 d h p r + k 32 d h p r O 3 d h p r + k 52 d h p r C 5 d h p r ( k 21 d h p r + k 23 d h p r ) O 2 d h p r ( k 24 d h p r + k 25 d h p r ) O 2 d h p r
(10)
d O 3 d h p r d t = k 23 d h p r O 2 d h p r k 32 d h p r O 3 d h p r
(11)
d C 4 d h p r d t = k 54 d h p r C 5 d h p r + k 24 d h p r O 2 d h p r ( k 45 d h p r + k 42 d h p r ) C 4 d h p r
(12)
d C 6 d h p r d t = k 36 d h p r O 3 d h p r k 63 d h p r C 6 d h p r
(13)
The open state O 3 dhpr accounts for the increased tail current produced as the result of a large depolarization. Expressions for rate constants:
k 42 d h p r = 1000.0
(14)
k 45 d h p r = 600000.0
(15)
k 54 d h p r = k 45 d h p r k 52 d h p r k 24 d h p r k 42 d h p r k 25 d h p r
(16)
k 15 d h p r = k 51 d h p r k 12 d h p r k 25 d h p r k 52 d h p r k 21 d h p r
(17)
Case 1: With CICR (control)
k 12 d h p r = ( 300.3808 301.0817 1 + e ( ν 13.5918 8.85 ) )
(18)
k 21 d h p r = ( 3359.8754 + 966.95 1 + e ( ν 1.66 1.4585 ) )
(19)
ξ = 550 + 6 × C a M K I I a c t + C a N a c t
(20)
k 24 d h p r = ( 336160 × S 2 d h p r ξ )
(21)
k 25 d h p r = ( 5939.4 + 306806.8 1 + e ( ν + 1.8716 1.3072 ) ξ )
(22)
k 52 d h p r = ( 0.02925 + 0.48961 1 + e ( ν + 12.2249 0.974 ) )
(23)
k 23 d h p r = ( 1352 1350 1 + e ( U 1.8 ) + 1700 1 + e ( v + 5 ) )
(24)
k 32 d h p r = ( 1030.0575 713.1966 1 + e ( v 5.0 ) )
(25)
k 51 d h p r = 383.5435 1 + e ( ν 8.3702 7.047 ) + 0.16725 1 + e ( ν 31.8252 0.01075 ) + 378.7085 1 + e ( ν 7.8953 6.7691 ) 373.015
(26)
k 36 d h p r = ( 100 100.0 1 + e ( ν 41.0 1.0 ) )
(27)
k 63 d h p r = ( 600.0 1 + e ( ν 41.0 1.0 ) )
(28)
Case 2: With Ryanodine applied (no RyR release)
k 12 d h p r = ( 30392.87 30394.73 1 + e ( ν 183.5975 28.5662 ) )
(29)
k 21 d h p r = ( 3568.74658 + 21921.25344 1 + e ( ν + 24.9838 0.7372 ) )
(30)
ξ = 550 + 6 × C a M K I I a c t + C a N a c t
(31)
k 24 d h p r = ( 336160 × S 2 d h p r ξ )
(32)
k 25 d h p r = 0.0893 + 26.3268 1 + e ( ν + 15.0 1.0 ) ξ
(33)
+ 32.6642 1 + e ( ν 14.4897 0.3131 ) + 63.808 1 + e ( ν 27.411 0.0854 ) ξ
(34)
k 52 d h p r = 2.69 + 0.74152 1 + e ( ν + 27.945 1.46976 ) 2.43098 1 + e ( ν 40.25 0.4034 )
(35)
k 23 d h p r = 1908 2687.3087 1 + e ( ν + 15.0 0.9998 ) 879.543 1 + e ( ν 11.1635 0.3238 )
(36)
k 32 d h p r = 145107.50934 150299.048 1 + e ( ν 1037.5224 312.91482 ) + 5537.7875 1 + e ( ν 9.4649 0.2893 )
(37)
k 51 d h p r = 27.651 1 + e ( ν 59.48 4.444 ) + 11.488 1 + e ( ν + 15.1148 0.9423 ) 16.3777 1 + e ( ν + 15.01 0.6351 ) + 5.3384
(38)
k 36 d h p r = ( 100 100.0 1 + e ( ν 41.0 1.0 ) )
(39)
k 63 d h p r = ( 600.0 1 + e ( ν 41.0 1.0 ) )
(40)
Case 3: With Ca2+ substituted with Ba2+
k 12 d h p r = 3028.7604 6279.833 1 + e ( ν 58.1484 9.5348 ) + 3259.706 1 + e ( ν 40.2856 0.2272 )
(41)
k 21 d h p r = 26296.0826 1 + e ( ν + 17.40398 0.22797 ) + 4718.8152 1 + e ( ν + 1.4892 4.9581 ) + 1544.9745 1 + e ( ν 132.7305 3.6234 ) 808.0319
(42)
k 24 d h p r = ( 0.0 )
(43)
k 25 d h p r = 5883.3476 + 5973.3014 1 + e ( ν 41.4221 0.5703 ) 528.171 1 + e ( ν 20.4274 0.187 )
(44)
k 52 d h p r = 175.78742 174.7872 1 + e ( ν 32.3072 4.4783 ) 229.8479 1 + e ( ν 36.9508 5.4521 )
(45)
k 23 d h p r = 2601.9597 1 + e ( ν 19.882 0.4459 ) 2647.5114 1 + e ( ν + 15.6872 2.1947 ) 2501.6394 1 + e ( ν 39.8404 0.3926 ) + 2647.2042
(46)
k 32 d h p r = 1683.7686 1 + e ( ν 20.0452 0.40326 ) + 2161.9699 1 + e ( ν 360.2862 1879.441 ) 5926.5225 1 + e ( ν + 1.7842 13.43183 ) + 4129.8882
(47)
k 51 d h p r = 11.1993 1 + e ( ν 35.2478 0.2466 ) + 36.2341 1 + e ( ν + 23.9615 15.8653 ) 12.163 1 + e ( ν 27.1529 0.9455 ) + 37.7673
(48)
k 36 d h p r = ( 100 100.0 1 + e ( ν 41.0 1.0 ) )
(49)
k 63 d h p r = ( 600.0 1 + e ( ν 41.0 1.0 ) )
(50)

Uptake of Ca2+ from the cytosol into the LSR

Ca2+ fluxes from cytosol to SERCA and SERCA to LSR
J c y t , s e r c a = k c y t , s e r c a × [ C a 2 + ] m y o 2 × S E R C A t o t k c y t , s e r c a × [ C a 2 + ] m y o 2 × [ C a 2 + ] s e r c a k s e r c a , c y t × [ C a 2 + ] s e r c a
(51)
J s e r c a , s r = k s r , s e r c a × [ C a 2 + ] L S R 2 × [ C a 2 + ] s e r c a k s r , s e r c a × [ C a 2 + ] L S R 2 × S E R C A t o t + k c y t , s r × [ C a 2 + ] m y o
(52)
I c y t , s e r c a = J c y t , s e r c a × 2 F V s e r c a
(53)
I s e r c a , s r = J s e r c a , s r × 2 F V s e r c a
(54)
Differential equation for Ca2+ buffered by the SERCA protein:
d [ C a 2 + ] s e r c a d t = J c y t , s e r c a J s e r c a , s r
(55)
Expressions for the rate constants for Ca2+ binding to/release from SERCA:
k c y t , s e r c a = k c y t , s e r c a ( 1.07 + 5500 × C a M K I I a c t )
(56)
K c y t , s e r c a = K c y t , s e r c a ( E C 50 f w d ) 2
(57)
k s e r c a , s r = K s e r c a , s r ( 1.07 + 5500 × C a M K I I a c t )
(58)
k s r , s e r c a = K s e r c a , s r ( E C 50 b w d ) 2
(59)
Affinities for the forward and backward Ca2+ fluxes:
E C 50 f w d = 0.015 ( 1 + P S R × P L B d p ) × ( 1.027 1 + 5500 × C a M K I I a c t )
(60)
E C 50 b w d = 1250 1110 × P S R × P L B d p
(61)
Differential equation for phospholamban
d P L B d p d t = k 12 P L B P L B p k 21 P L B ( 5500 C a M K I I a c t + 2.5 C a N a c t + P K A a c t ) 2 P L B d p
(62)
P L B p = 1 P L B d p
(63)

Ca2+ pump in SL

I P M C A = I P M C A ¯ ( [ C a 2 + ] m y o k m p c a + [ C a 2 + ] m y o )
(64)

Na+/Ca2+ exchanger

I N a C a = N U M 1 × ( N U M 2 N U M 3 ) ( 1 + 0.27 e 0.65 F V R T )
(65)
N U M 1 = R N a C a ( V max 1 + ( K m A l l o [ C a 2 + ] N a C a ) 2 )
(66)
N U M 2 = e 0.35 F V R T ( [ N a + ] N a C a ) 3 [ C a 2 + ] o
(67)
N U M 3 = e 0.65 F V R T ( [ N a + ] o ) 3 [ C a 2 + ] N a C a
(68)
= K m C a o ( [ N a + ] N a C a ) 3 + K m N a o 3 [ C a 2 + ] N a C a + K m N a i 3 [ C a 2 + ] ( 1 + [ C a 2 + ] N a C a K m C a i ) + K m C a i ( [ N a + ] o ) 3 ( 1 + ( [ N a + ] N a C a K m N a i ) 3 ) + ( [ N a + ] N a C a ) 3 [ C a 2 + ] o + ( [ N a + ] o ) 3 [ C a 2 + ] N a C a
(69)

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 .

Na+/Cs+ pump

I N a C s = R N a C s I N a C s ¯ ( [ C s + ] 0 [ C s + ] 0 + k m c s ) × ( 100 ( [ N a + ] m y o ) 1.5 ( [ N a + ] m y o ) 1.5 + k m n a 1.5 ) × ( 0.65558 0.18445 + e ( ν + 53.353 15.58 ) )
(70)

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

I N a , b = G N a b ( ν R T F ln ( [ N a + ] o [ N a + ] m y o ) )
(71)

Ca2+ transfer from LSR to a single jSR

i t r = ( [ C a 2 + ] L S R [ C a 2 + ] j S R τ t r ) 2 F V j S R
(72)

Ca2+ release from a unit jSR into a single DCU

i r y r = J r y r × ( 2 F × π r 2 z )
(73)
where,
J r y r = O 2 r y r ( [ C a 2 + ] j S R [ C a 2 + ] r y r ) P r y r π r 2 z
(74)
Gating scheme for the 4-state Markovian model for the RyR-sensitive SR Ca2+ release channel:
C 4 r y r = 1 ( C 1 r y r + O 2 r y r + C 3 r y r )
(75)
d C 1 r y r d t = k 41 r y r C 4 r y r + k 21 r y r O 2 r y r ( k 12 r y r + k 14 r y r ) C 1 r y r
(76)
d O 2 r y r d t = k 12 r y r C 1 r y r + k 32 r y r C 3 r y r ( k 21 r y r + k 23 r y r ) O 2 r y r
(77)
d C 3 r y r d t = k 23 r y r O 2 r y r + k 43 r y r C 4 r y r ( k 34 r y r + k 32 r y r ) C 3 r y r
(78)
k 12 r y r = [ C a 2 + ] r y r 2 ( 0.05 3.7465 × 10 2 1 + e [ C a 2 + ] r y r 2052.7 5.983 ) v a r + ( C a M K I I a c t 12000 )
(79)
k 21 r y r = ( 698.56 618.56 1 + e [ c a 2 + ] r y r 2052.9 4.35 ) × ( 1 + 6.0 [ C a 2 + ] r y r ) 1 v a r
(80)
k 14 r y r = 8.748 × 10 2 [ C a 2 + ] r y r v a r
(81)
k 41 r y r = ( 2.0 v a r ) + C a M K I I a c t 12000
(82)
k 43 r y r = k 12 r y r ; k 34 r y r = k 21 r y r ;
(83)
k 23 r y r = k 14 r y r ; k 32 r y r = k 41 r y r ;
(84)
v a r = [ 10.0 × e A 1 l s 0.7013 0.03 ] 2
(85)
where [Ca2+] ryr is the Ca2+ concentration at the mouth of the RyR channel on the dyadic side. In the presence of caffeine (CF, concentration in μ M):
O 2 r y r = ( 0.522 1 + e [ C F ] 410 264.9 ) + 0.503
(86)
Gating scheme for the 6-state Markovian model for the Luminal Calcium sensor.
I 4 l s = 1 ( A 1 l s + I 2 l s + I 3 l s )
(87)
B 5 l s = 1 ( I 2 l s + I 3 l s + B 6 l s )
(88)
d A 1 l s d t = k 41 l s I 4 l s + k 21 l s I 2 l s ( k 12 l s B 5 l s + k 14 l s ) A 1 l s
(89)
d I 2 l s d t = k 12 l s B 5 l s A 1 l s + k 42 l s B 5 l s I 4 l s + k 32 l s I 3 l s ( k 21 l s + k 24 l s + k 23 l s + k 25 l s ) I 2 l s
(90)
d I 3 l s d t = k 43 l s B 5 l s I 4 l s + k 23 l s I 2 l s ( k 34 l s + k 32 l s ) I 3 l s
(91)
d B 6 l s d t = k 56 l s [ C a 2 + ] j S R B 5 l s k 65 l s B 6 l s
(92)

A 1 ls    Fractional occupancy of RyR by Triadin/Junctin

I 2 ls    Fractional occupancy of RyR by Triadin/Junctin and Calsequestrin

I 3 ls    Fractional occupancy of Triadin/Junctin by Calsequestrin

I 4 ls    Free Triadin/Junctin, B5ls - Free Calsequestrin

B 6 ls    Fractional Occupancy of Calsequestrin by Calcium

Ca2+ dependent CaM mediated activation of CaMKII and CaN

Ca2+ binding to CaM and CaM buffering
C a M = C a M t o t C a 2 C a M C a 4 C a M C a M B C a 2 C a M B C a 4 C a M B C a M C a N C a 2 C a M C a N C a 4 C a M C a N C a M K I I t o t ( P 1 + P 3 + P 5 + P 4 )
(93)
B = B t o t ( C a M B + C a 2 C a M B + C a 4 C a M B )
(94)
R O 2 = k 02 C M [ C a 2 + ] 2 × C a M k 02 C M × C a 2 C a M
(95)
R 24 = k 24 C M [ C a 2 + ] 2 × C a 2 C a M k 42 C M × C a 4 C a M
(96)
R O 2 B = k 02 B C M [ C a 2 + ] 2 × C a M B k 20 B C M × C a 2 C a M B
(97)
R 24 B = k 24 B C M [ C a 2 + ] 2 × C a 2 C a M B k 42 B C M × C a 4 C a M B
(98)
R O B = k 0 B o n C M × C a M × B k 0 B o f f C M × C a 2 C a M B
(99)
R 2 B = k 2 B o n C M × C a 2 C a M × B k 2 B o f f C M × C a 2 C a M B
(100)
R 4 B = k 4 B o n C M × C a 4 C a M × B k 4 B o f f C M × C a 4 C a M B
(101)
d C a 2 C a M d t = R 02 R 24 R 2 B R 2 C a N + C a M K I I t o t × ( R C K 56 R C K 21 )
(102)
d C a 4 C a M d t = R 24 R 4 C a N R 4 B C a M K I I t o t × ( R C K 46 R C K 23 )
(103)
d C a 4 C a M d t = R 24 R 4 C a N R 4 B C a M K I I t o t × ( R C K 46 R C K 23 )
(104)
d C a M B d t = R 0 B R 02 B
(105)
d C a 2 C a M B d t = R 02 B + R 2 B R 24 B
(106)
d C a 4 C a M B d t = R 24 B R 4 B
(107)
CaMKII activation:
P 2 = 1 P 3 P 1 P 4 P 5 P 6
(108)
T = P 3 + P 4 + P 5 + P 6
(109)
k 34 C K = 0.055 × T + 0.0074 × T 2 + 0.015 × T 3
(110)
R C K 34 = k 34 C K × P 3 ( K P P 1 × P P 1 t o t × P 4 k m p p 1 + ( C a M K I I t o t × P 4 ) )
(111)
R C K 21 = k 21 C K × C a 2 C a M × P 2 k 12 C K × P 1
(112)
R C K 13 = k 13 C K × [ C a 2 + ] 2 × P 1 k 32 C K × P 3
(113)
R C K 23 = k 23 C K × C a 4 C a M × P 2 k 32 C K × P 3
(114)
R C K 45 = k 45 C K × P 4 k 54 C K × [ C a 2 + ] 2 × P 5
(115)
R C K 46 = k 46 C K × P 4 k 64 C K × C a 4 C a M × P 6
(116)
R C K 56 = k 56 C K × P 5 k 65 C K × C a 2 C a M × P 6
(117)
R C K 51 = k P P 1 × P P 1 t o t × P 5 k m P P 1 + ( C a M K I I t o t × P 5 )
(118)
R C K 62 = k P P 1 × P P 1 t o t × P 6 k m P P 1 + ( C a M K I I t o t × P 6 )
(119)
d P 1 d t = R C K 21 + R C K 51 R C K 13
(120)
d P 3 d t = R C K 23 + R C K 13 R C K 34
(121)
d P 4 d t = R C K 34 R C K 46 R C K 45
(122)
d P 5 d t = R C K 45 R C K 56 R C K 51
(123)
d P 6 d t = R C K 46 R C K 56 R C K 62
(124)
C a M K I I a c t = 100 × ( P 3 + P 4 + P 5 + P 6 )
(125)
CaN activation:
C a 2 C a N = C a N t o t C a 4 C a N C a M C a N C a 2 C a M C a N C a 4 C a M C a N
(126)
R C N C a 4 = k C a o n C N × [ C a 2 + ] 2 × C a 2 C a N k C a o f f C N × C a 4 C a N
(127)
R C N 02 = k 02 C N × [ C a 2 + ] 2 × C a M C a N k 20 C N × C a 2 C a M C a N
(128)
R C N 24 = k 24 C N × [ C a 2 + ] 2 × C a 2 C a M C a N k 42 C N × C a 4 C a M C a N
(129)
R C N 0 = k 0 o n C N × C a M × C a 4 C a N k 0 o f f C N × C a M C a N
(130)
R C N 2 = k 2 o n C N × C a 2 C a M × C a 4 C a N k 2 o f f C N × C a 2 C a M C a N
(131)
R C N 4 = k 4 o n C N × C a 4 C a M × C a 4 C a N k 4 o f f C N × C a 4 C a M C a N
(132)
d C a 4 C a N d t = R C N C a 4 R C N 0 R C N 2 R C N 4
(133)
d C a M C a N d t = R C N 0 R C N 02
(134)
d C a 2 C a M C a N d t = R C N 2 + R C N 02 R C N 24
(135)
d C a 4 C a M C a N d t = R C N 4 + R C N 24
(136)
C a N a c t = 100 ( C a 4 C a M C a N + 0.1 C a 2 C a M C a N ) + 10 C a 4 C a M C a N + 0.1 C a M C a N +0 .1 C a 4 C a N
(137)

A2 - Differential equations for buffers used in the model

Fluorescent indicator dye
d [ C a F 3 ] d t = k f l u o 3 + [ C a m y o 2 + ] ( [ f l u o 3 ] t o t [ C a F 3 ] ) k f l u o 3 [ C a F 3 ]
(138)

Intracellular Ca2+ buffering:

calmodulin (bulkmyoplasm):
d O c d t = 200 [ C a 2 + ] m y o ( 1 O c ) 476.0 × O c
(139)

Troponin

(Fractional occupancy of troponin-Ca complex by Ca2+):
d O t c d t = 78.4 [ C a 2 + ] m y o ( 1 O t c ) 392.0 × O t c
(140)
(Fractional occupancy of troponin-Mg complex by Ca2+):
d O c d t = 200 [ C a 2 + ] m y o ( 1 O t m g c O t m g m g ) 6.6 × O t m g c
(141)
(Fractional occupancy of troponin-Mg complex by Mg2+):
d O t m g m g d t = 2.0 [ M g 2 + ] m y o ( 1 O t m g c O t m g m g ) 666.0 × O t m g m g
(142)

A3 - Differential equations for ion concentrations used in the model

Intracellular Ca 2+ concentration:

1. Ca 2+ concentration in the cytosol
d [ C a 2 + ] m y o d t = I d y a d I P M C A I C a , L , S L 2 F V m y o + I c y t , s e r c a + 2 I N a C a , S L 2 F V m y o d O d t d C a F 1 d t
(143)
where I dyad is the net integrated Ca2+ flux diffusing out of all the dyadic units into the cytosol
d O d t = 3.2 d O T C d t + 6.4 d O T M g C d t + 1.8 d O c d t
(144)
2. Ca 2+ concentration in the jSR
d [ C a 2 + ] j S R d t = i t r i r y r 2 F V j S R 31000 d B 6 l s d t
(145)
3. Ca 2+ concentration in the LSR
d [ C a 2 + ] L S R d t = I c y t , s e r c a N d y a d i t r 2 F V L S R
(146)
4. Diffusion equation for Calcium in the dyadic space
[ C a 2 + ] d y a d t = D C a 2 [ C a 2 + ] d y a d + J r y r + J d h p r + J N a C a + J b n d
(147)
J d h p r = ( i C a , L , T T N d y a d ) ( 5.1821 × 10 3 π r 2 z )