This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.106.089425v1 92/10/3379    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Tanskanen, A. J. Articles by Winslow, R. L. Search for Related Content PubMed PubMed Citation Articles by Tanskanen, A. J. Articles by Winslow, R. L.

Protein Geometry and Placement in the Cardiac Dyad Influence Macroscopic Properties of Calcium-Induced Calcium Release

Antti J. Tanskanen ^*  , Joseph L. Greenstein ^*  , Alex Chen ^*  , Sean X. Sun   and Raimond L. Winslow ^*  

^* The Institute for Computational Medicine, Center for Cardiovascular Bioinformatics and Modeling, The Whitaker Biomedical Engineering Institute, Department of Biomedical Engineering, and Department of Mechanical Engineering, The Johns Hopkins University School of Medicine and Whiting School of Engineering, Baltimore, Maryland

Correspondence: Address reprint requests to Joseph L. Greenstein, Clark Hall, Rm. 201D, 3400 N. Charles St., Baltimore, MD 21218. E-mail: jgreenst{at}jhu.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
In cardiac ventricular myocytes, events crucial to excitation-contraction coupling take place in spatially restricted microdomains known as dyads. The movement and dynamics of calcium (Ca²⁺) ions in the dyad have often been described by assigning continuously valued Ca²⁺ concentrations to one or more dyadic compartments. However, even at its peak, the estimated number of free Ca²⁺ ions present in a single dyad is small (10–100 ions). This in turn suggests that modeling dyadic calcium dynamics using laws of mass action may be inappropriate. In this study, we develop a model of stochastic molecular signaling between L-type Ca²⁺ channels (LCCs) and ryanodine receptors (RyR2s) that describes: a), known features of dyad geometry, including the space-filling properties of key dyadic proteins; and b), movement of individual Ca²⁺ ions within the dyad, as driven by electrodiffusion. The model enables investigation of how local Ca²⁺ signaling is influenced by dyad structure, including the configuration of key proteins within the dyad, the location of Ca²⁺ binding sites, and membrane surface charges. Using this model, we demonstrate that LCC-RyR2 signaling is influenced by both the stochastic dynamics of Ca²⁺ ions in the dyad as well as the shape and relative positioning of dyad proteins. Results suggest the hypothesis that the relative placement and shape of the RyR2 proteins helps to "funnel" Ca²⁺ ions to RyR2 binding sites, thus increasing excitation-contraction coupling gain.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Contraction of the cardiac myocyte results from a process known as excitation-contraction coupling (ECC). ECC is initiated when individual L-type calcium (Ca²⁺) channels (LCCs) open in response to membrane depolarization, producing Ca²⁺ flux into a microdomain known as the dyad. The resulting increase in dyadic Ca²⁺ leads to opening of Ca²⁺-sensitive Ca²⁺-release channels known as ryanodine receptors (RyR2s) located in the closely apposed junctional sarcoplasmic reticulum (JSR) membrane, producing additional flux of Ca²⁺ from the JSR into the dyad (Fig. 1 A). These two sources of Ca²⁺ flux generate the intracellular Ca²⁺ transient that triggers cardiac muscle contraction. Understanding the molecular basis of this Ca²⁺-induced Ca²⁺-release (CICR) process is therefore of fundamental importance to understanding cardiac muscle function in both health and disease.

View larger version (33K): [in this window] [in a new window]   FIGURE 1  (A) A schematic representation of calcium-induced calcium release in the dyad. Ca2+ ions pass through L-type Ca2+ channels (LCCs) and enter the dyad where they bind to ryanodine receptors (RyRs), triggering their opening leading to Ca2+ release from the junctional sarcoplasmic reticulum (JSR). Ca2+ ions released from the SR diffuse out of the dyad and into the bulk myoplasm, where they can bind to myofilaments, and initiate cell shortening. (B) Layout of LCCs and RyR2s within the dyad. The quasicrystal array of RyR2s (gray shaded boxes with crosshair) lie in the JSR membrane with a smaller number of LCCs (green circles) randomly distributed in the opposing T-tubule membrane.

The small dimensions, in particular the limited height, of the dyad imply that the structure of proteins located within the dyad may also serve to restrict motion of Ca²⁺ ions. The largest protein within the dyad is RyR2. The structure of RyR2 has been measured at a resolution of 3.0 nm (7), whereas that of RyR1 (the skeletal muscle isoform, sharing 70% sequence identity to RyR2 (8)) has been measured at 1.0 nm resolution (9). RyR2 is a large protein composed of four 565-kDa subunits. The cytoplasmic portion of the protein has dimensions 27 x 27 x 12 nm (9,10), where the 12-nm height spans nearly the full distance between JSR and T-tubule membranes (1,2,11). The clustering of RyR2s opposite each LCC therefore presents a considerable Ca²⁺ diffusion barrier. The crystal structure of the cardiac LCC has also been measured at 0.4-nm resolution (12). The LCC is 19 nm in height and 14.5 nm in width, protruding from the T-tubule membrane 2 nm into the dyad (12). In addition, the structure of the Ca²⁺-binding protein calmodulin (CaM), one molecule of which is tethered to the inner pore of the LCC (13), has been measured at 0.1-nm resolution (14). CaM is a dual-lobed protein of approximate dimensions 4.5 x 4.5 x 6.5 nm. Given the small dimensions of the dyad, it is likely that the physical location and dimensions of these important dyad proteins have a considerable influence on motion of Ca²⁺ ions and thus properties of CICR.

In this study, we present a computational, molecularly and structurally detailed model of the cardiac dyad and use this model to address the following questions: 1), What is the number of Ca²⁺ ions that are present in the dyad during CICR? 2), How does the physical arrangement of large proteins within the dyad influence the process of CICR? 3), How does "signaling noise" due to the small number of Ca²⁺ ions within the dyad affect the nature of CICR. The results demonstrate that: 1), at the level of the single dyad, CICR may be mediated by as few as 20–50 Ca²⁺ ions; 2), even though the number of free Ca²⁺ ions in a single dyad during CICR is small and highly variable, ECC gain is consistently high (measured over the population of dyads in a single cell) and is a robust feature of local Ca²⁺ signaling; and 3), the structure of proteins that reside in the dyad help to funnel Ca²⁺ toward RyR2 binding sites, and in doing so, enhance ECC gain.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
The motion of individual Ca²⁺ ions in the dyad is influenced by the following factors: a), interactions between Ca²⁺ ions and other mobile ions and molecules within the dyad; b), the physical boundary imposed by membranes bounding the dyad and the location and shape of proteins within the dyad; c), the presence of an electric field near the T-tubule membrane; d), stochastic gating of LCCs and RyR2s, producing a stochastic boundary flux; e), location of Ca²⁺ binding sites (such as those on LCCs, RyR2s, and CaM); and f), the nature of the interface between the dyad and the surrounding myoplasm. To properly simulate Ca²⁺ dynamics, these influences must be included in a detailed model of the dyad.

Structure of the dyad
The volume of the dyad is represented using a 200 x 200 x 15 nm lattice (Fig. 1 B; gray boxes represent RyR2; green dots represent LCCs) with 1-nm spacing between lattice points. The lower boundary of the dyad is limited by the T-tubular membrane and the upper boundary is limited by the JSR membrane (Fig. 1 A). Ca²⁺ can diffuse across the lateral boundary of the dyad between the dyadic volume and the myoplasm (Fig. 1 A). The geometry of the dyadic volume that is accessible to diffusing Ca²⁺ ions is determined by the presence, shape, and location of proteins, which due to their large size relative to the dyad, occupy a significant fraction of the dyadic volume. Chief among these are LCCs, RyR2s, and CaM. RyR2s are located in the JSR membrane in quasicrystalline arrays of tens of RyR2s (Fig. 1 B) (1). In this study, we employ a structural model of the cytosolic (i.e., dyadic) assembly of the RyR2 based on experimental cryo-electron microscopy (cryo-EM) measurements (Fig. 2 A) (10,15). The overall dimensions of the cytosolic assembly of the RyR2 are 27 x 27 x 12 nm and the model exhibits fourfold symmetry with four "feet" structures. Substantial evidence indicates the pore of the RyR2 complex is located in the center of the tetrameric structure (Fig. 2 A, green dot) (9). Each model dyad is assumed to contain 20 RyR2s arranged in an asymmetric 4 x 5 quasicrystal (Fig. 1 B), where neighboring RyR2s come into contact with each other near their corners, with an overlap of 12 nm along their edges (16). A 5:1 RyR2/LCC ratio is assumed consistent with previous measures of the relative density of RyR2 and LCC proteins (2,17). Thus, this model dyad shares some fundamental features with the calcium release unit model studied previously by Greenstein and Winslow (18).

View larger version (32K): [in this window] [in a new window]   FIGURE 2  (A) Model RyR2 structure as seen from its cytosolic face based on experimentally measured structure (10,15,28). Red spheres indicate the positions of the baseline model Ca2+-binding activation sites, blue spheres indicate the positions of alternate hypothetical activation sites, yellow spheres indicate the positions of model inactivation sites, and the green sphere depicts the assumed position of the RyR2 pore. (B) Side view of a portion of the dyad showing a single RyR2 and opposing LCC with tethered CaM. The assumed position of the Ca2+-binding sites for CDI on CaM are indicated by yellow spheres.

View larger version (38K): [in this window] [in a new window]   FIGURE 6  Histograms of Ca2+ release latency (defined as the time between opening of the first LCC and opening of the first RyR2) from a population of 1500 dyads simulated at 0 mV under the following conditions: (A) geometric protein structures included, electric field present (i.e., membrane surface charges present); (B) geometric protein structures excluded, electric field present; (C) geometric protein structures included, electric field absent; and (D) geometric protein structures excluded, electric field absent.

The locations of LCCs within the dyad are shown in Fig. 1 B (green dots). The cytosolic region of each LCC structure occupies an area of 10 x 13 nm and extends 2 nm from the inner surface of the T-tubule membrane. A constitutively tethered CaM acts as the Ca²⁺ sensor for Ca²⁺-dependent inactivation of LCCs (31,32). It has recently been shown that a single CaM molecule is both necessary and sufficient to produce Ca²⁺-dependent inactivation of its associated channel (13). The model therefore includes a single CaM molecule associated with each LCC (Fig. 2 B). Because the precise location and orientation of this CaM with respect to the LCC is not known, we assume, for simplicity, that the CaM is located adjacent to the LCC along the sarcolemma. A geometric model of CaM is created by approximating the crystal structure of its Ca²⁺-unbound form (identified as 1CFD in Protein Data Bank), as measured experimentally by Kuboniwa et al. (33), at 1-nm resolution (Fig. 2 B). A single CaM molecule contains four Ca²⁺ binding sites, of which the two sites located on the carboxyl tail have been shown to be responsible for Ca²⁺-dependent inactivation (CDI) of the associated LCC (31,32). In this model, CDI can proceed when both of the carboxyl-tail binding sites are occupied (Fig. 2 B, yellow dots). It has been assumed that the carboxyl tail of CaM and the associated CDI binding sites are oriented toward the LCC. LCC facilitation is thought to be regulated by the two Ca²⁺ binding sites located on the amino-terminal lobe of CaM (34) (see Anderson (35) for a review), however the process of LCC facilitation was not included in this model, and the amino-terminal binding sites were therefore not implemented in the model CaM molecule.

RyR2 and LCC kinetic models
A comprehensive description of Ca²⁺ dynamics in the dyad requires kinetic models of LCC and RyR2 gating to quantitatively describe the source fluxes of Ca²⁺ into the dyad. LCCs and RyR2s are modeled using continuous-time, discrete-state Markov processes. In many previous models (e.g. (18,36)), the incorporation of Ca²⁺-dependent ion channel gating kinetics was accomplished by allowing model state transition rates to be defined by expressions that depend upon the relevant Ca²⁺ concentration. This approach is an approximation in that it combines the Ca²⁺-binding step and the conformational change of the ion channel (e.g., inactivation) into a single state transition where the conformational change of the channel protein has been assumed to be rate-limiting (e.g., see Peterson et al. (32)). In the model presented here, Ca²⁺ binding to the channel protein must be considered separately from the subsequent conformational change of the channel protein. In order for a Ca²⁺-dependent state transition to occur, the required Ca²⁺-binding sites must first be occupied by Ca²⁺ ions (e.g., Ca²⁺-dependent inactivation of an LCC requires that the two carboxyl-terminal Ca²⁺ binding sites of the associated CaM are occupied). When a freely diffusing Ca²⁺ ion enters a lattice position adjacent to an available Ca²⁺-binding site, that Ca²⁺ ion has the opportunity to bind to the site. The relative magnitude of the binding rate compared to the rate of diffusion determines the probability that the Ca²⁺ ion either binds to the site or diffuses away. Ca²⁺-binding transitions are incorporated into the overall model of Ca²⁺ movement in the dyad (see below). Ion channel transition rates are therefore defined as functions that depend upon the occupancy of the Ca²⁺ binding sites (rather than Ca²⁺ concentration). Rates for binding and unbinding of Ca²⁺ to sites on both the LCCs and RyR2s are given in the Appendix.

The kinetic model of the LCC is described using an 11-state continuous time Markov chain model (Fig. 3 A) developed previously (18,36,37). Briefly, the upper row of states describes the LCC normal gating mode (Mode Normal) and the lower row of states describes gating when the LCC undergoes Ca²⁺-dependent inactivation (Mode CDI). The "downward" transition rates from Mode Normal to Mode CDI become nonzero only when both Ca²⁺-binding sites of the associated CaM molecule are occupied. The rate for CDI was adjusted based on that used in a previous implementation of this channel model in the presence of 100 µM [Ca²⁺] (18). Transition rates from Mode CDI to Mode Normal (i.e., recovery from Ca²⁺-mediated inactivation) do not depend on whether or not Ca²⁺ ions are bound to the associated CaM. Voltage-dependent inactivation of the LCC is incorporated as a separate gating variable (not shown in Fig. 3) with rates identical to those described previously (18).

View larger version (7K): [in this window] [in a new window]   FIGURE 3  (A) Kinetic model of LCC gating based on the model of Greenstein and Winslow (18). (B) Kinetic model of RyR2 gating based on that of Stern et al. (38). The processes of Ca2+ binding/unbinding to activation and inactivation sites and the process of voltage-dependent inactivation are not depicted in this figure (see text for details). Transition rate constants are provided in the Appendix.

The unitary Ca²⁺ current of a single RyR2 is assumed to be 1.24 pA under physiological conditions (38,39). This corresponds to an influx rate of 3870 Ca²⁺ ions per millisecond. LCC permeability is set such that the unitary Ca²⁺ current through a single LCC is 0.12 pA at 0 mV (41), which corresponds to an influx rate of 375 Ca²⁺ ions per millisecond. LCC open-channel permeability, and hence Ca²⁺ influx rate, varies with membrane potential as described previously (see Fig. 2 G of Greenstein and Winslow (18)). The entry of Ca²⁺ ions into the dyad via an open RyR2 or LCC is simulated by a Poisson process with a rate set equal to the rate of Ca²⁺ ion entry for each LCC or RyR2 as described above. It is assumed that ion entry via an LCC or RyR2 is blocked if a Ca²⁺ ion occupies the lattice point that is adjacent to the channel mouth.

Membrane Ca²⁺ binding
Ca²⁺ buffering plays an important role in cardiac myocyte Ca²⁺ dynamics (2). In the dyad, membrane phospholipid headgroups act as fixed Ca²⁺ buffers. The density of slow SR and sarcolemmal buffering sites is large, typically assumed to be 0.1–0.2 site/nm² (3,4). Binding site density and Ca²⁺ binding/unbinding rates in this model are based on the work of Langer and Peskoff (3) and Soeller and Cannell (4) and are given in the Appendix. It has been assumed that both low-affinity and high-affinity Ca²⁺ binding sites are present on both the SR membrane and sarcolemma. Ca²⁺ binding to buffer sites is implemented in the same manner as described above for Ca²⁺ binding to ion channel proteins. In addition to fixed Ca²⁺ buffers, mobile Ca²⁺ buffers are also known to be present in the dyad. Freely diffusing calmodulin is the main mobile Ca²⁺ buffer in the dyad, however, the typically assumed CaM concentration of 15 µM (42) translates into a single molecule in a dyad of radius 50 nm. Recently, however, it has been suggested that there may be local enrichment of CaM molecules in the vicinity of the LCC, resulting in as many as 25 free CaM molecules at the site of each LCC (13). It remains unclear how these CaM molecules would be targeted to the LCC, and whether they can freely diffuse in this restricted space. Because the local dynamics of free CaM are not yet understood, and because CaM is a large, rather slowly diffusing molecule, freely diffusing CaM molecules were not included in this dyad model. Similarly, the role of ATP as a mobile Ca²⁺ buffer was not included in this model.

Ca²⁺ dynamics in the dyad
The motion of a mobile Ca²⁺ ion in the dyad is influenced by the Brownian random force from the surrounding solvent and the electrostatic potential stemming from proteins, membranes, and other ions, including other Ca²⁺ ions. Indeed, in the confined space of the dyad, collisions of Ca²⁺ ions with other free ions are likely to be frequent and therefore correlations between the ions should be considered. Consequently, we model the joint positions (r₁,...,r_N) of N Ca²⁺ ions present in the dyad as a 3N-dimensional Brownian motion in a potential field, which describes both interactions of Ca²⁺ ions with other Ca²⁺ ions as well as electrostatic potentials. The time evolution of the joint probability density of these Ca²⁺ ions to be present at positions (r₁,...,r_N) in the dyad at time t, P(r₁,...,r_N,t), is described by the Fokker-Planck equation (FPE; see, e.g., Risken (43))

(1)

where r_i = (r_i,1,r_i,2,r_i,3) is a vector indicating the position of the i-th ion, D is the diffusion constant, k_B is Boltzmann's constant, T is temperature, and the notation /r_i is defined as

The total potential energy of the system, V, is given as a function of the ion positions,

(2)

where q is the elementary charge. The total potential energy V has several contributions: a), (r_i) is the electrostatic potential of the i-th ion due to charges on the surrounding lipids and proteins (for example, contains the potential due to surface charge density, , from the negatively charged phospholipid headgroups); b), u(r_i) is known as a hard-core potential that becomes nonzero at the location of impenetrable structures such as proteins within the dyad. This potential simply defines the space accessible for the mobile ions and is ultimately determined by the structural model of the dyad; and c), U(r₁,...,r_N) is the mutual interaction potential between the Ca²⁺ ions. This potential is determined by the physical size (hard-core repulsion) of the ions, and the dielectric/buffering conditions in the dyad. For two ions, U depends on the ionic separation and the range of U is determined by the Debye length, , within the dyad. For Ca²⁺ ions in the dyad, is 1 nm (4). In effect, the ions feel a strong repulsion if they are within distance and no interaction otherwise. Therefore serves as a natural correlation length between the ions.

The electrostatic potential, , is dominated by membrane surface charges in the dyad. We use the Debye-Hückel model of charge-charge interaction, in which is given by

(3)

where the integral is taken over the surface S of the membrane, is the membrane surface charge density, and is again the Debye length. The constant ₀ depends on the dielectric constant and ionic conditions in the dyad. We follow Soeller and Cannel (4) and approximate as a monoexponential function arising from sarcolemmal surface charges

(4)

where h is the height of the dyad and z is the vertical distance from the sarcolemma at point r (see Fig. 1 of Soeller and Cannell (4)).

The volume within the dyad that is accessible to a particular Ca²⁺ ion is assumed to include any space not occupied by an LCC, RyR2, CaM, or other Ca²⁺ ions. It is assumed that diffusing Ca²⁺ ions cannot penetrate the surface of any protein (including the RyR2 foot processes), the JSR membrane, or the sarcolemma. These surfaces/membranes are therefore considered as reflecting (i.e., no-flux boundary conditions). The lateral boundary at the interface between the dyadic volume and myoplasm is treated as an absorbing boundary for ions flowing from the dyad to the myoplasm (see below). In addition, a baseline rate of Ca²⁺ entry into the dyad from the myoplasm has been implemented based on the assumption of a constant myoplasmic Ca²⁺ concentration of 100 nM. These boundary conditions ensure that the dyads are statistically independent.

Discrete approximation
The multidimensional FPE (Eq. 1) describes the evolution of the probability density function for the positions of Brownian particles subject to a potential. Rather than solving the time-dependent joint probability density from the FPE, we generate paths of Ca²⁺ ion movement in the dyad using an algorithm described by Wang et al. (44). This method produces a finite differencing of the FPE that can be interpreted as a spatially discrete Markov process, which can then be simulated using Monte Carlo methods (45). The dyadic space is discretized into a lattice consisting of 1 nm³ boxes. Because the timescale of Ca²⁺ diffusion is sufficiently more rapid than other kinetic events such as channel gating, the system can be considered to be in a local steady state within any single box. Furthermore, if the potential changes linearly within the box, then an analytic local solution is possible within the box (44,46). Using the local solution, and continuity requirements, Brownian motion in continuous space can be discretized into a discrete-state Markov process. Consequently, the FPE is converted into a master equation

(5)

where = (r₁, r₂,..., r_N) labels the composite state of the system where each r_i describes the position of the i-th ion. P is the probability that the system is in state , and is the transition rate from state to state '. The transition rates for Ca²⁺ movement can be directly obtained from the local equilibrium solution of the FPE (44). Transitions can only occur between connected states of the system (i.e., states that differ in the position of a single ion by a distance equivalent to a single lattice point). Therefore if = (r₁, r₂,...,r_i,..., r_N) and ' = (r₁, r₂,...,r_i',..., r_N), then

(6)

where is the transition rate from the initial location r_i of the i^th ion to the final location r_i' of the i^th ion. Transition rate formulas for a variety of boundary conditions have been derived previously (44) as:

(7)

where D is the diffusion constant of Ca²⁺, is the lattice spacing (= 1 nm), and is the change in the local potential given by

(8)

Note that if there is no potential difference between two adjacent boxes, , which is equal to the inverse of the expected time for the ion to diffuse to the adjacent box in the absence of an electric field. Under conditions where would be large and negative (e.g., when an ion occupies the adjacent box or there is an adjacent reflecting boundary), is automatically set to zero. Thus, no more than one Ca²⁺ ion can occupy a single lattice point at a time and Ca²⁺ ions cannot transition into inaccessible region of the dyad (e.g., space filled by protein structures). On the lateral boundaries of the simulation grid, absorbing boundary conditions are applied (Eq. 7), allowing ions to escape the dyadic region and be absorbed into the bulk cytosolic compartment. These definitions completely specify the transformation from the continuous-space FPE description to the discrete-space Markov description.

In discretized form, the description of Ca²⁺ dynamics in the dyad can be coupled with the processes of Ca²⁺ binding and gating dynamics of RyR2s and LCCs. As described above, gating dynamics for each channel are described by a Markov process where the binding and unbinding of Ca²⁺ ions may affect ion channel transition rates. The Markov processes for Ca²⁺ diffusion, Ca²⁺ binding, and channel gating can be combined. Thus, the state space described by is expanded to include both the positions of Ca²⁺ ions and the states of the ion channels. The process of Ca²⁺ binding to the channels is incorporated by allowing ions in lattice boxes adjacent to a binding site to bind with rates k_on and k_off (see Appendix for rates governing each type of binding site). The incorporation of these processes allows the entire dyad to be simulated by a single all-inclusive Markov process, in which the dynamics of the RyR2s and LCCs are coupled to the dyadic Ca²⁺ dynamics of the model.

To compute the time-dependent solution, we perform a kinetic Monte Carlo sampling algorithm to generate paths of mobile ions. The general algorithm of sampling continuous time Markov processes was first proposed by Bortz et al. (45). This method has been used extensively in the study of biochemical networks and molecular motors (44,46,47). Briefly, for each given state , all possible destination states, ', and their associated transition rates, are tabulated. The net exit rate for the current state is obtained by summing over ',

(9)

A random number t is then picked from the waiting time distribution and is assigned as the time the system will remain in its current state. To determine the destination state for the following transition, another uniform random number in the interval (0,K) is generated. The subinterval in which resides, where the size of each subinterval corresponds to the magnitude of each individual exit rate from the current state, determines the destination state. In this fashion, trajectories/paths that include ion diffusion, opening and closing of channels, and binding of Ca²⁺ to channels can be computed. Time-dependent distributions of Ca²⁺ ions can be obtained by analyzing sample paths generated using the Monte Carlo procedure. Average quantities such as ion current, Ca²⁺ flux, and ECC gain are computed by summing/averaging over a large number of dyad simulations. In addition to calculating average quantities, the Monte Carlo approach used here allows for analysis of the degree of variance in relevant measures such as ECC gain.

A typical simulation involves an ensemble of 400 dyads, each of which is clamped to test potentials in the range of –40 mV to +50 mV in 10-mV increments where each voltage clamp step is held for 100 ms. This simulation requires 40 h of runtime on an IBM BladeCenter with 62 nodes, each with two Intel Xeon 2.8 Ghz processors. Although the computational task of simulating this model is substantial compared to most conventional myocyte models, the runtime is several orders of magnitude faster than a typical molecular dynamics simulation. This model allows for simulation of a large number of dyads with timescales measured in seconds. A 64-bit version of the Mersenne Twister algorithm (48,49) was employed as the random number generator.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX ACKNOWLEDGEMENTS REFERENCES

 
Excitation-contraction coupling gain
Macroscopic SR Ca²⁺ release is often quantified by measuring ECC gain. The characteristic voltage-dependent shape of the ECC gain function arises as a result of local control of Ca²⁺ release at the level of the dyad (18). ECC gain measures the relative magnitude of JSR Ca²⁺ release via RyR2s to Ca²⁺ influx via LCCs. Fig. 4 A shows the voltage dependence of peak LCC Ca²⁺ flux (J_LCC, solid circles) and peak RyR2 Ca²⁺ flux (J_RyR, open circles) obtained in response to 100-ms duration depolarizing voltage-clamp steps from –40 to 50 mV from a holding potential of –100 mV. In Fig. 4 B, the peak fluxes of Fig. 4 A have been normalized based on their respective maxima. Although both peak Ca²⁺ fluxes are essentially bell shaped, J_RyR reaches its peak value at a potential that is 10 mV more negative than that where J_LCC reaches its peak (i.e., the curves are shifted with respect to each other), as was also observed in experiments (50,51). The consequence of this shift is that peak ECC gain, defined as the ratio of peak J_RyR to peak J_LCC, is a monotonically decreasing function of membrane potential, as shown in Fig. 4 C. In Fig. 4 C, the simulated peak ECC gain is compared to the experimental results of Song et al. (50) and Wier et al. (51). The simulated peak ECC gain curve exhibits gain values that are within the range of the experimentally measured values and demonstrates a similar shape. Only at potentials more negative than –20 mV do the simulation results differ from the experimental measurements. However, this is not unexpected because both experimental and model (see below) measures of ECC gain exhibit higher variability in this range of potentials compared to positive potentials. In Fig. 4 D, integrated ECC gain (defined here as the ratio of the total number of Ca²⁺ ions that enter the dyad via RyR2s to the total number that enter via LCCs over the duration of a 100-ms voltage clamp) is shown for comparison to peak ECC gain. Both measures of gain display the characteristic monotonic decay with increasing membrane potential. Each data point in the panels of Fig. 4 was calculated as the average of three simulations, each consisting of a population of 400 dyads.

View larger version (13K): [in this window] [in a new window]   FIGURE 4  Voltage dependence of LCC Ca2+ influx (JLCC), JSR Ca2+ release flux (JRyR), and ECC gain. (A) Peak values of JLCC (solid circles) and JRyR (open circles) as a function of clamp membrane potential. (B) Peak Ca2+ fluxes (data of panel A) normalized by their respective maxima. (C) ECC gain based on peak Ca2+ fluxes (data of panel A) for the model (solid line) are compared to experimental data of Song et al. (dark gray dashed line (50)) and Wier et al. (light gray dashed line (51)). (D) Comparison of peak ECC gain (solid line) and integrated ECC gain (dashed line).

View larger version (14K): [in this window] [in a new window]   FIGURE 5  The role of dyadic protein structures on ECC gain. Peak ECC gain as a function of membrane potential for the baseline model, which includes space-filling geometric models of protein structure in the dyad (solid line), for the model with protein structures excluded (dashed line), and for a modified model with dyad height reduced from 15 to 13 nm and protein structures excluded (dotted line). Three independent simulations of gain including (circles) and excluding (triangles) geometric protein structures are shown at –20, 0, and +20 mV.

The above results indicate that protein structures are an important factor that influence the dynamics of CICR, and suggest that the location of Ca²⁺-binding sites on these structures might also affect CICR. To explore this possibility, an alternate set of Ca²⁺ activation binding sites on the RyR2 were tested in the model. On each subunit, the alternate activation site was located on the surface of domain 1, within 3 nm of the RyR2 pore (Fig. 2 A, blue dots). Fig. 7 shows that with the alternate activation sites (dashed line, average of three simulations of 400 dyads each), peak ECC gain is increased nearly twofold compared to the baseline model (solid line). The increase in gain follows from the fact that the LCCs are aligned directly across the dyadic cleft from RyR2s such that each LCC pore is located directly across from a RyR2 pore. Hence, the alternative activation binding sites are located significantly closer to the source of Ca²⁺ trigger flux than in the control case, increasing the probability that Ca²⁺ ions encounter and bind activation sites on the RyR2. In a manner similar to that demonstrated earlier for the control RyR2 binding site locations, the RyR2 protein structure plays an important role in funneling the Ca²⁺ ions to the alternative Ca²⁺-binding sites. The exclusion of protein structures still leads to decreased ECC gain in the case where RyR2 Ca²⁺-binding sites have been moved to the alternative locations (e.g., with alternative binding site locations, at 0 mV ECC gain decreases to 14.5 from 19.7 upon removal of protein structures). This finding demonstrates that various locations along the surface of the RyR2s encounter different numbers of Ca²⁺ ions, and this suggests that RyR2 activity, and hence efficiency of CICR, may vary considerably with changes in the alignment LCCs and RyR2s.

View larger version (13K): [in this window] [in a new window]   FIGURE 7  The role of RyR2 Ca2+-binding sites and sarcolemmal surface charges on peak ECC gain. Peak ECC gain as a function of membrane potential for the baseline model (solid line), for the model with Ca2+-binding activation sites relocated near the RyR2 pore (dashed line), and for the model with no membrane surface charges (dash-dotted line).

Number of Ca²⁺ ions in dyad
The peak number of free Ca²⁺ ions in the dyad has been estimated to be small, 1 free Ca²⁺ ion at a Ca²⁺ concentration of 10 µM in a dyad of radius 50 nm (2). Fig. 8 demonstrates the contribution of LCCs and RyR2s to the number of free Ca²⁺ ions in the local vicinity of each channel type within a dyad. The results of Fig. 8 A are obtained from a representative simulation of a model dyad containing only a single LCC in the sarcolemma located at the center of the dyad, where the number of Ca²⁺ ions within a 50-nm radius of the LCC is shown. The membrane potential is clamped to 0 mV, the initial state of the LCC is open (gray bars at top indicate when the LCC is open), and the initial number of Ca²⁺ ions is zero. The number of Ca²⁺ ions is sampled every 0.1 ms (i.e., at 10 kHz) for a duration of 10 ms. In this dyad, the average number of Ca²⁺ ions within 50 nm of the LCC is 0.46 (dashed gray line), however, there is a high degree of variability in the number present at any given instant (mean ± SD 1.2, dotted gray line). In this example, the maximum number of Ca²⁺ ions present is seven, whereas at most times there are zero Ca²⁺ ions in the dyad. Because each Ca²⁺ ion corresponds to a Ca²⁺ concentration of 14.1 µmol/L (within the 50-nm radius of the LCC), the equivalent Ca²⁺ concentration is effectively changing in large discrete steps as Ca²⁺ ions enter and exit this volume, reaching a maximum of 98.7 µmol/L corresponding to seven Ca²⁺ ions. Fig. 8 B shows the average number of Ca²⁺ ions sampled in the vicinity of an LCC calculated over 1000 independent simulations using the same protocol. At the beginning of the protocol, there is an average of 2 Ca²⁺ ions in the vicinity of the LCC. However, as a result of CDI, LCC open probability decreases and the average number of Ca²⁺ ions drops to 0.5 after 3 ms. These results suggest that the number of free Ca²⁺ ions in the dyad arising from LCC activity is sufficiently small and variable such that their description as a continuously varying Ca²⁺ concentration may not be adequate.

View larger version (27K): [in this window] [in a new window]   FIGURE 8  Counting the number of free Ca2+ ions in the dyad arising from a single LCC or RyR2. In each panel, the number of Ca2+ ions is sampled every 0.1 ms from the volume within a cylinder of 50 nm radius centered at the single LCC or RyR2. (A) Representative simulation of a dyad containing only a single gating LCC (initially open) with membrane potential clamped at 0 mV, showing the number of free Ca2+ ions (black bars), average (dashed gray line), and standard deviation (dotted gray line). The gray bars at the top indicate when the LCC is open. (B) Mean number of free Ca2+ ions based on 1000 dyad simulations of the protocol described in panel A. (C) Representative simulation of a dyad containing only a single gating RyR2 (initially open), showing the number of free Ca2+ ions (black bars), average (dashed gray line), and standard deviation (dotted gray line). The gray bars at the top indicate when the RyR2 is open. (D) Mean number of free Ca2+ ions based on 1000 dyad simulations of the protocol describe in panel C.

CICR in the dyad
Panels A and B of Fig. 9 demonstrate fundamental features of a representative Ca²⁺ release event in the dyad stimulated via a voltage clamp to 0 mV at time zero. Fig. 9 A shows the number of free Ca²⁺ ions in the dyad resulting from LCC and RyR2 gating, and Fig. 9 B shows the number of LCCs (gray line) and RyR2s (black line) that are open as a function of time during this release event. In this example, two LCCs open at 3 ms following the voltage clamp and the Ca²⁺ influx into the dyad triggers the opening of three RyR2s 1–2 ms later, one of which inactivates after 7 ms, the next inactivates after an additional 2 ms, and the final one inactivates 13 ms after the start of the release event. The temporal shape of the Ca²⁺ signal in the dyad (Fig. 9 A) is closely correlated to the number of open RyR2s (Fig. 9 B). Fig. 9 C demonstrates an average dyadic Ca²⁺ signal, as calculated by averaging the number of Ca²⁺ ions at each point in time for 1000 independent dyad simulations. As would be expected the peak of the average signal is reduced compared to that of an individual dyad due to temporal dispersion of Ca²⁺ release events and the fact that not all dyads will exhibit a Ca²⁺ release event. The duration of the average local Ca²⁺ transient shown in Fig. 9 C is 20 ms (at half-maximal amplitude), similar to that measured for Ca²⁺ spikes in myocytes using confocal imaging techniques (50,53,54). Fig. 9 D demonstrates the snapshot of the spatial distribution of Ca²⁺ ions in a single dyad (as a function of the distance from the center) during a release event similar to that shown in Fig. 9 A. The Ca²⁺ ion density was calculated as a function of radial distance from the dyad center as the mean density taken over a 5-ms time window centered at 1.4 ms following the beginning of the release event (a time at which two RyR2s were open). In this case the open RyR2s are centrally located in the dyad, and the density of Ca²⁺ ions is highest at this location and decreases with increasing radius from the release site.

View larger version (34K): [in this window] [in a new window]   FIGURE 9  A CICR event in a single dyad. The response is evoked by a 0-mV voltage clamp step occurring a time 0. (A) The number of free Ca2+ ions in the dyad volume as a function of time. (B) The number of open RyR2s (black solid line) and the number of open LCCs (gray solid line) in the dyad during the release event depicted in panel A. (C) The average number of free Ca2+ ions in the dyad as calculated from 1000 independent dyads subject to the same voltage clamp protocol as the event shown in panel A. (D) Average density (ions nm–3) of Ca2+ ions as a function of radial distance from the center of the dyad at 1.4 ms following a release event. (E) Distribution of the duration of Ca2+ release events. (F) Distribution of the maximum number of RyR2s that open simultaneously during a single Ca2+ release event (bars, left axis) and the mean (± SD) of the maximal number of open LCCs (taken over the 5-ms window preceding each release event) associated with each grouping of release events (solid circles, right axis).

Ca²⁺-binding sites are present in the sarcolemmal and JSR membranes of the dyad and have been included in this model as described in Methods. Fig. 10 A shows the ratio of free Ca²⁺ ions to total Ca²⁺ ions, calculated as the average ratio over 317 simulated release events occurring in response to a voltage clamp to 0 mV at time zero. This ratio is initially 17% (on average), corresponding to the early phase of Ca²⁺ release. Immediately following release, Ca²⁺ ions diffuse away from the LCC and RyR2 structures and many encounter and bind to available Ca²⁺ buffering sites. As a result, the ratio of free Ca²⁺ ions decreases to an equilibrium value of 2%, as can be seen at times >10 ms following Ca²⁺ release in Fig. 10 A. For the same set of release events, Fig. 10 B shows the average fraction of free Ca²⁺ ions that are located within 1 nm of the sarcolemmal membrane. Within 1 ms of the beginning of the clamp to 0 mV, and for the entire duration of the clamp, 29–30% of the free calcium ions tend to be restricted to the space adjacent to the sarcolemmal membrane as a result of the electric field. This is consistent with the finding that ECC gain is increased in the absence of the electric field (Fig. 7), as the electric field reduces the likelihood that free Ca²⁺ ions will encounter RyR2 binding sites. The accumulation of Ca²⁺ ions adjacent to the sarcolemma is qualitatively consistent with the steep gradient in Ca²⁺ concentration reported adjacent to the sarcolemma in the continuum model of Soeller and Cannell (4).

View larger version (9K): [in this window] [in a new window]   FIGURE 10  The effects of Ca2+ buffers and membrane surface charges on free Ca2+ in the dyad. (A) The average ratio of the number of free Ca2+ ions to the total number of Ca2+ ions during CICR is shown as a function of time. (B) The average fraction of free Ca2+ ions located within 1 nm of the sarcolemmal membrane. The responses in both panels are evoked by a 0-mV voltage clamp step occurring at time 0, and the average shown is obtained from 317 release events.

View larger version (29K): [in this window] [in a new window]   FIGURE 11  Variability of ECC gain. (A) Peak ECC gain as a function of voltage is shown for 10 repeated model simulations based on a population of 400 dyads. (B) Mean peak ECC gain (black solid line) based on the data of panel A is shown along with the mean ± SD (denoted S.D., dashed lines). (C) The distribution of peak ECC gain at –40 mV for the whole myocyte (12,500 dyads) computed using bootstrap replicates. (D) The distribution of peak ECC gain at 0 mV for the whole myocyte computed using bootstrap replicates.