This Article Abstract Full Text (PDF) Supplemental All Versions of this Article: biophysj.105.065466v1 90/6/1999    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 Koh, X. Articles by Levchenko, A. Search for Related Content PubMed PubMed Citation Articles by Koh, X. Articles by Levchenko, A.

A 3D Monte Carlo Analysis of the Role of Dyadic Space Geometry in Spark Generation

Whitaker Institute for Biomedical Engineering and Department of Biomedical Engineering, The Johns Hopkins University Whiting School of Engineering, Baltimore, Maryland

Correspondence: Address reprint requests to Andre Levchenko, Whitaker Institute for Biomedical Engineering, Dept. of Biomedical Engineering, The Johns Hopkins University Whiting School of Engineering, 208C Clark Hall, 3400 North Charles St., Baltimore, MD 21218. Tel.: 410-516-5584; Fax: 410-516-6240; E-mail: alev{at}jhu.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
In multiple biological systems, vital intracellular signaling processes occur locally in minute periplasmic subspaces often referred to as signaling microdomains. The number of signaling molecules in these microdomains is small enough to render the notion of continuous concentration changes invalid, such that signaling events are better described using stochastic rather than deterministic methods. Of particular interest is the dyadic cleft in the cardiac myocyte, where short-lived, local increases in intracellular Ca²⁺ known as Ca²⁺ sparks regulate excitation-contraction coupling. The geometry of dyadic spaces can alter in disease and development and display significant interspecies variability. We created and studied a 3D Monte Carlo model of the dyadic cleft, specifying the spatial localization of L-type Ca²⁺ channels and ryanodine receptors. Our analysis revealed how reaction specificity and efficiency are regulated by microdomain geometry as well as the physical separation of signaling molecules into functional complexes. The spark amplitude and rise time were found to be highly dependent on the concentration of activated channels per dyadic cleft and on the intermembrane separation, but not very sensitive to other cleft dimensions. The role of L-type Ca²⁺ channel and ryanodine receptor phosphorylation was also examined. We anticipate that this modeling approach may be applied to other systems (e.g., neuronal growth cones and chemotactic cells) to create a general description of stochastic events in Ca²⁺ signaling.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
In the past decade, it has become apparent that multiple intracellular signaling processes occur in minute periplasmic subspaces, often referred to as signaling microdomains (1,2). These are membrane-restricted subcellular spaces where vital, local signaling events can take place to regulate cell function. In these small reaction volumes (in the attoliter range) (3–7), the number of reactant molecules corresponding to physiologically relevant concentrations can be exceedingly low. The reaction system thus cannot be accurately described by a continuous, deterministic approach based on reaction-diffusion "partial differential equations", since fluctuations around average values become important and stochastic behavior dominates (8). Here, the familiar macroscopic notion of concentration is no longer useful, necessitating the use of stochastic methods to better describe signaling events. One important example of a microdomain where stochastic events are expected to be prominent is the dyadic cleft of the ventricular myocyte. The dyadic cleft spans an estimated 10- to 12-nm gap (3–7,9) between the voltage-gated L-type Ca²⁺ channels (LCCs) on the transverse tubule (TT) membrane, and ryanodine receptors (RyRs) on the sarcoplasmic reticulum (SR). In response to depolarization of the cardiac action potential (AP), Ca²⁺ influx via LCCs traverses the intracellular compartment to trigger Ca²⁺ release from the closely apposed SR via adjacent Ca²⁺-sensitive RyRs (Fig. 1 A). These events can trigger a fundamentally important positive feedback mechanism of amplifying Ca²⁺ signals known as Ca²⁺-induced Ca²⁺ release (CICR) (10).

View larger version (47K): [in this window] [in a new window]   FIGURE 1  Visualization of the model layout and elements. (A) Extracellular Ca2+ (red) enters the dyadic space via LCCs (blue) to trigger SR Ca2+ release from the Ca2+-sensitive RyRs (yellow). (B) The standard model cleft is assumed to be a rectangular space where the distance between the TT and SR membranes, denoted as Dcleft, is normally 12 nm. The lateral section is square with the x,y dimensions denoted as Lcleft. Bearing in mind the potential variability of cleft geometry, we compare models with fixed numbers of channels and Lcleft varying from 100 to 500 nm. Depicted in the figure is the case where Lcleft = 100 nm. (C) Here, the channel density remains unchanged, whereas the model dimensions are changed. The dimensions are Lcleft,1 = 100 nm, Lcleft,2 = 400 nm, and Dcleft = 12 nm. (D) Models of varying Dcleft are compared: 12 nm (left), 25 nm (center), and 40 nm (right). Lcleft is fixed at 100 nm. An increase in the average distance between LCCs and RyRs may reduce the triggering of SR Ca2+ release.

Over the years, there has been mounting evidence that the control of CICR is contingent upon local Ca²⁺ in the immediate vicinity of the channels, rather than whole-cell Ca²⁺ (12–14). In 1992, Stern proposed the local control theory of EC coupling (see, e.g., Greenstein and Winslow (4), Stern (13), and Niggli and Lederer (15)), which asserts that tight regulation of CICR is made possible by the clustering of LCCs and RyRs into discrete units of Ca²⁺ release, rendering them sensitive to local rather than global Ca²⁺ levels. In other words, macroscopic Ca²⁺ release events are intrinsically controlled by the conductance properties (gating and ion permeation) of individual LCCs and RyRs, and the relative spatial localization of the two channel types. The trigger LCC Ca²⁺ influx has been found to have a tight, smoothly graded control of SR Ca²⁺ release, despite the fact that CICR is a self-regenerative process that intuitively leads to an all-or-none response. This paradox of Ca²⁺ regulation can be explained by the spatial coupling of LCCs and RyRs into functional release units that are able to operate independently, hence providing further evidence for the local control theory. In the rest of this article, we will refer to the site of local control, or the functional Ca²⁺ release unit, as the dyadic cleft (4).

In addition to the type and number of Ca²⁺ channels in the dyadic cleft, the geometry of this functionally significant microdomain is critical to the analysis of Ca²⁺ regulation. The dyadic cleft has been estimated to have a radius of 0.05–0.2 µm, and a height of 10–12 nm (3–7), which gives rise to a volume 1.0–20.0 x 10^–13 µL if a cylindrical cleft geometry is assumed. During Ca²⁺ spark formation, the local [Ca²⁺] in this small compartment may reach a peak of 100 µM to 1 mM during the first 10 ms, and decline rapidly with a halftime of 20 ms to the diastolic [Ca²⁺] of 1 µM or less (3,6,11,16). These concentration values translate to a maximum of 1000 ions at peak [Ca²⁺], and only 1–2 Ca²⁺ ions at rest. Hence, it is clear that at the level of local Ca²⁺ signaling, the number of Ca²⁺ ions can be small enough to render the notion of continuous concentration changes invalid (8). One aspect of cleft geometry that may be a possible feature in failing myocytes is a dramatic phenotypical change in structure, which has been associated with abnormal Ca²⁺ signaling (17–19).

Another important aspect of Ca²⁺ regulation is the possibility that various signaling events can be capable of modulating the properties of molecular species responsible for Ca²⁺ release and uptake. In microdomains like the dyadic cleft, the ß-adrenergic receptor (ß-AR) signaling complex may be colocalized with immediate downstream effectors into multimolecular complexes for specific signal transduction to Ca²⁺ channels (20,21). Numerous studies have shown a correlation between sustained ß-adrenergic stimulation and cardiac dysfunction (22). It is also known that ß-adrenergic agonists are able to increase force development and accelerate contractile relaxation by altering channel gating properties. Excessive ß-adrenergic stimulation, however, may lead to LCC and RyR hyperphosphorylation, which can impair contractility and cause heart failure (22–24).

In recent years, several computational models incorporating either fully deterministic or partially deterministic methods (3,4,7,25–27) have been developed to explore Ca²⁺ dynamics in the dyadic space. However, Ca²⁺ signaling at the level of individual clefts is highly stochastic and dependent on local Ca²⁺, suggesting that nondeterministic simulation methods are necessary to achieve a sufficiently detailed description of signaling events.

Here we use a 3D Monte Carlo approach to better describe Ca²⁺ regulation within dyadic clefts of cardiac myocytes. To account for events that involve individual diffusing and reacting molecules, we examined a completely stochastic simulation of localized Ca²⁺ signaling. Simulations are performed with MCell (www.mcell.cnl.salk.edu), a software package designed for realistic simulations of cellular microphysiology. Using this approach, we are able to track the positions of individual molecules on different spatial and temporal scales, as well as predict how specific changes in the physical environment might influence Ca²⁺ signaling. We incorporate, to the largest possible extent, information from available literature into a dyadic cleft model of hypothetical dimensions. Despite extensive studies over the years on Ca²⁺ spark properties, the exact nature and definition of sparks still remain controversial (28,29). Our model is intended to capture the properties of the dyadic cleft as a fundamental site of Ca²⁺ release (9,30–32) by describing in detail the basic elements thought to be necessary for spark formation. In addition, we consider the role of geometry and protein kinase A (PKA)-mediated phosphorylation in Ca²⁺ regulation in the cleft.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Structural considerations
The modeled cleft space
Imaging studies in skeletal and cardiac muscles have shown that clusters of RyRs are physically separated on junctional domains of the SR membrane (jSR) (9). It has also been reported that, depending on species, a few tens to up to 200 RyRs are clustered in regular lattice arrays (7). A recent study on channel locations in frog myocytes (33) showed that LCCs also reside in clusters on the sarcolemmal membrane. Moreover, optical experiments examining the distribution of proteins involved in EC coupling indicate that LCCs and RyRs are colocalized into discrete Ca²⁺ release units (CaRU) (34,35). Dissipation of local [Ca²⁺] from each CaRU takes place so rapidly that it does not normally trigger SR Ca²⁺ release at neighboring junctions. The above findings provide strong evidence for the spatial and functional independence of single dyadic clefts, and support the notion that discretized units of Ca²⁺ release microdomains play a major role in EC coupling. In keeping with these observations, we simulate individual dyadic clefts separately as independent units.

Cleft geometry
For the purposes of our model we assume that the standard cleft is a rectangular space with a square lateral section of x,y dimensions denoted as L_cleft, and a z dimension of 12 nm between the TT and SR membranes, denoted as D_cleft (Fig. 1 B). Ca²⁺ that escapes from the cleft is unlikely to return due to rapid dissipation into the larger surrounding cytosolic space. Hence the cleft sidewalls are regarded as Ca²⁺-absorbing surfaces, whereas the top and bottom walls are Ca²⁺-reflecting to account for the impermeability of cell membranes (more on boundary conditions below).

The separation between the SR and exterior membranes may be relatively conserved at 10–12 nm; however, other geometrical properties of clefts differ widely across varying species and muscle types (9). In previous cardiac models, the dyadic cleft has been estimated to span 0.05–0.2 µm in radius and 10–12 nm in height (3–7). However, the complex microarchitecture, more specifically the shape and size of this space, is not explicitly known. Bearing in mind the potential variability of cleft geometry, we compare models with L_cleft ranging from 0.1–0.5 µm, and volumes ranging from 1.2–30.0 x 10^–13 µL. We compare different channel densities by fixing the numbers of channels, and also compare clefts of different dimensions where the channel density remains unchanged (Fig. 1 C).

Boundary conditions
Although the dyadic cleft is modeled as a functionally independent unit, it is not isolated in the sense that we do take into account the escape of Ca²⁺ into the surrounding cytosol. MCell maps the positions of surfaces and effector sites (e.g., channel proteins and receptors) in space, and tracks the positions of individual diffusing ligands at every time step. When a ligand is detected at a point of intersection with a surface, there can be different possible outcomes depending on the properties of the surface. For example, at the point of intersection, the surface may be reflective, transparent, absorptive, or occupied by an effector site with an associated chemical reaction mechanism. After entering the dyadic cleft, free Ca²⁺ ions may escape into the neighboring space by diffusion. To account for this escape, the four sidewalls of the rectangular cleft are modified to absorb any Ca²⁺ ions that come into contact. This modification is based on the assumption that the return of Ca²⁺ into the cleft by diffusion is insignificant, since the cleft volume is considerably smaller than the surrounding cytosolic volume. The top and bottom surfaces are modified to reflect Ca²⁺ back into the cleft, to represent the impermeable phospholipid bilayers of the TT and SR membranes, respectively.

The placement of absorptive surfaces at the boundaries of the cleft raised the possibility of a discontinuity in the diffusion gradient along the edges of the cleft. To investigate this effect, large volumes were placed around the cleft, which moved the absorptive surfaces away from the boundaries of the cleft and allowed reentry of Ca²⁺ into the cleft. It was found that the introduction of a pericleft volume to the model did not significantly affect the characteristics of the original Ca²⁺ response (see Supplementary Material). This affirms that the return of Ca²⁺ into the cleft is insignificant and that absorptive boundary conditions can be assumed.

Disruption of TT morphology
There is evidence that TT morphology transforms during myocyte development and heart failure, and that these transformations might be closely linked to observed functional variations in the cell. Instances of such morphology changes in pathologic conditions include abnormally shaped TTs in hypertrophic human heart (17), and TT damage or loss in canine tachycardia-dilated cardiomyopathy (18).

TTs are composed of interconnected elements resembling caveolae (19), which are infoldings of the surface membrane in cells lacking TTs (e.g., atrial cells and neonatal cardiomyocytes). In a study on Ca²⁺ sparks produced in caveolae (36), it was proposed that a 20- to 50-nm increase in the distance between the membranes on which LCCs and RyRs reside might produce sparks with altered spatiotemporal characteristics.

When Wang et al. applied a G-seal patch-clamp on the TT membrane, the efficacy of CICR was observed to be greatly reduced, most likely due to a significant degree of LCC/RyR uncoupling caused by membrane deformation associated with the G-seal (37). This anomalous phenomenon was not observed when a loose M-seal patch-clamp was applied instead. Studies involving aligned and nonaligned myocytes have also associated changes in myocyte morphology with altered channel properties (38), possibly due to the disrupted spatial relationship between channels.

The observations described above point toward the likelihood that disruption of the dyadic cleft geometry, including D_cleft, can interfere with the tight local control between LCCs and RyRs, and lead to pathological alterations in Ca²⁺ signaling. A spatial change in LCC/RyR organization in failing hearts may adversely affect the triggering of SR Ca²⁺ release by LCC Ca²⁺ influx (39). This spatial change is often associated with an increase in the average distance between LCCs and RyRs (40). Although the structure of our model does not closely replicate the 3D microarchitecture of a dyadic cleft, we are able to take advantage of its simplicity to make predictions about the effects of altering D_cleft. We first observe the Ca²⁺ release events that take place in a model with a normal D_cleft of 12 nm. We then observe the effects of changing D_cleft systematically to small (9 nm) or large (up to 40 nm) values existing outside the physiological range believed to be possible (Fig. 1 D).

Model formulation
Channel stoichiometry
Due to the numerous unknowns in the system we are modeling, it is challenging, if not impossible, to create a model that can reproduce physiological events with absolute accuracy. Next we describe the experimental data serving as the basis for some of the assumptions we make.

When modeling the dyadic cleft, one should keep in mind that the reported densities of Ca²⁺ channel proteins are often ambiguous since they are dependent on the intrinsic architecture of the clefts and may also vary between different species (41) and muscle types. In an ultrastructural study on rat ventricular myocytes (42), the densities of LCCs and RyRs are reported to be 84/µm² and 765/µm², respectively. In the body muscles of arthropods, the density of RyRs has been estimated to be 1275–1890/µm² (41), which translates to 15–28 RyRs residing on a square patch of jSR membrane area of L_cleft 0.1–0.15 µm (CaRUs in the body muscles of arthropods and the skeletal and cardiac muscles of vertebrates have been found to have similar architectures (41)). In specifying the number of channels in each cleft, we consider the above-mentioned data on channel densities, as well as the likely stoichiometry of LCC/RyR coupling. It was reported that a single LCC is able to trigger the opening of four to six RyRs to generate a spark (37), using a patch-clamp method in conjunction with confocal microscopy. In most simulations described below, a hypothetical dyadic cleft of variable L_cleft is composed of five LCCs and 25 RyRs.

LCC gating and permeation
There have been extensive studies in recent years that analyze and characterize single LCC properties in healthy heart cells, and also under a range of conditions that might indicate an increased susceptibility to disease (23,43). We model the stochastic gating of LCCs by specifying unimolecular transitions between adjacent open and closed states. MCell calculates reaction transition probabilities according to user-specified rate constants and determines which transitions will occur per time step. Rate constants (units s^–1) for LCC gating are derived from a study on single LCC availability and open probability (P_o) in healthy and failing human ventricles (43) (Table 1).

For each channel, only one reaction event (e.g., single ion entry and channel gating) may occur per Monte Carlo decision time step. Although the channel is in an open state, Ca²⁺ influx takes place via transitions of the open state back to itself (Fig. 2, A–D), according to a unimolecular transition probability that is calculated from the specified rate of ion generation. Unitary current can be modeled with reasonable accuracy if the simulation time step is sufficiently small. For example, 1.3 pA of current is equivalent to the entry of four divalent ions per µs, or one ion per 0.25 µs. This implies that a time step of 0.25 µs is small enough to simulate any current not exceeding 1.3 pA. For the purposes of our model, we use a 0.25-µs time step throughout all simulations. This time step has been shown to be sufficiently small, as demonstrated by additional studies in which a range of smaller time steps (up to 10 times smaller at 0.025 µs) were used (see Supplementary Material).

View larger version (21K): [in this window] [in a new window]   FIGURE 2  Complete state diagrams for the LCC (A) and RyR (C and D). These models incorporate a common schematic for channel gating and permeation that is summarized in panel B. The generation of Ca2+ (indicated by asterisks) is facilitated by open-state (O) transitions back to themselves according to a probability of unimolecular transition. (A) The LCC enters a Ca2+-inactivated state after two Ca2+-binding events associated with high-affinity CaM sites, whereas recovery from Ca2+ inactivation is governed by dissociation of Ca2+ from the CaM/LCC complex. Voltage inactivation may proceed from any of the LCC-residing states, and is assumed to be irreversible during the simulation of 250-ms duration. Both nonphosphorylated and phosphorylated LCC models follow the scheme depicted in panel A, but use different reaction parameters that were obtained from Schröder et al. (43) (see Table 1). (C) The nonphosphorylated RyR model is adopted from a Markovian model described by Saftenku et al. (49), where we also obtained the parameters associated with the closed states C1–C5, and the open states O1–O3. The implementation of RyR Ca2+-dependent inactivation is depicted by dotted lines, and proceeds upon Ca2+ binding to handpicked states (O1, C3) found to be prevalent throughout the simulations. (D) We introduce a new model for the phosphorylated RyR (see Methods), which is similar in principle to that depicted in Fig. 2 B, with the addition of a Ca2+-dependent inactivation process. Parameters are obtained from a study by Uehara et al. (23) on gating kinetics of phosphorylated RyR (see Table 2).

Inactivation of LCC currents
The inactivation of LCC currents serves as a way to constrain trigger Ca²⁺ influx during the cardiac AP and allows the EC coupling system to return to basal conditions during diastole. It is controlled by processes dependent on membrane potential and intracellular Ca²⁺ (46,47,44), known as voltage inactivation and Ca²⁺ inactivation, respectively. When the LCC is held persistently at a depolarized membrane potential, it can enter and stay in the voltage-inactivated state until membrane repolarization is introduced. In this model, the LCC is able to probabilistically enter a voltage-inactivated state from every reaction state (Fig. 2 A), governed by a rate constant derived from experimental findings (44) (see Table 1).

Ca²⁺ inactivation is a negative feedback mechanism triggered by elevated levels of intracellular Ca²⁺. It has been established that calmodulin (CaM), which is constitutively tethered to the LCC complex, is the Ca²⁺ sensor for this mode of inactivation (48). Results from the same study also suggest that Ca²⁺-dependent inactivation is contingent on the binding of two Ca²⁺ ions to each of the high-affinity sites of CaM, after which the Ca²⁺/CaM complex undergoes a conformational change that leads to channel inactivation. The LCC remains in this nonconducting state until the inactivation process is reversed in the event of a decrease in intracellular Ca²⁺ level. The Ca²⁺-binding and unbinding events leading to and from Ca²⁺ inactivation are reflected in the detailed model schematic (Fig. 2 A; (45)).

RyR state diagram
Saftenku et al. have assessed and described the gating behavior of single cardiac RyRs, using maximum likelihood analysis to estimate single-channel kinetic parameters from experimentally observed dwell times (49). Based on their studies (49), we adopt a Markovian model of the cardiac RyR that reportedly ranked highest among the others, according to the Schwartz criterion. The adopted schematic accounts for activation of the RyR by Ca²⁺, as well as channel gating between adjacent closed and open states. To complete the description of RyR activity, we append an additional reaction state that represents Ca²⁺-dependent inactivation (Fig. 2 C).

RyR permeation
We incorporate a unitary RyR current (i_RyR) of 0.1 pA, equivalent to an ion generation rate 0.35 x 10⁶ s^–1. As is the case for LCC unitary current, this value for i_RyR falls within the range of ion generation that may be accurately represented using a 0.25-µs time step. It has been reported that physiological unitary RyR current should be <0.6 pA (50), and possibly as little as 0.07 pA (16). The unitary current value we use corresponds to that from the formulation of an earlier cardiac model (4), where the concentration-dependent SR release flux through single RyRs can be approximated to give rise to unitary currents of 0.09–0.15 pA during a cardiac cycle.

Termination of SR Ca²⁺ release
Intuitively, CICR seems to be a self-regenerating process operating on the positive feedback of SR release Ca²⁺. This scenario implies the possibility of unstable global Ca²⁺ oscillations, which we know is not the case in healthy myocytes. A negative control mechanism must exist to terminate the Ca²⁺ spark by interrupting SR Ca²⁺ release. In 1985, Fabiato first proposed a Ca²⁺-dependent, negative feedback process of CICR inactivation (10). In later years, other mechanisms were proposed, and there has been a growing consensus that the negative control mechanism is actually a composite of processes acting in concert to terminate Ca²⁺ release through the RyRs (7,53). The different hypotheses include (we refer to reviews in (51,52)):

1. Ca²⁺-dependent inactivation. Experiments have shown that RyRs in cardiac cells are sensitive to jSR Ca²⁺, and turn off after free [Ca²⁺] reach high levels.
   
2. Stochastic attrition. All the RyRs in a release unit have the inherent quality of closing randomly. This process is sensitive to the number of channels in a cluster (13).
   
3. Adaptation.
   
4. Depletion of jSR or lumenal [Ca²⁺].
   
5. Time-dependent RyR inactivation.
   

The exact negative control mechanism is still a matter of debate. We choose to model Ca²⁺-dependent inactivation as a likely mechanism for the termination of SR Ca²⁺ release. Ca²⁺-dependent inactivation is incorporated into the RyR gating scheme by introducing state transitions into an inactivated state as mediated by Ca²⁺ binding (Table 2 and Fig. 2 C) (53). Such transitions are able to proceed from hand-picked states (O1, C3) found to be prevalent in the scheme. Once inactivated, the RyR is able to recover and return to the basal state via Ca²⁺ dissociation. In a parameter sensitivity study shown in Supplementary Material, we find the range of rates for RyR recovery from inactivation under which the system can display stable behavior. Should RyR recovery happen too readily, there is a possibility that premature sparks can occur before the onset of the next AP. This implies a deregulation of the EC coupling mechanism and a possible trigger of arrhythmic events. Although the addition of the Ca²⁺-inactivated state provides a pathway that was not present in the original model of Saftenku et al. (49), we find this modification necessary to appropriately describe pertinent properties, in particular the decline of elementary Ca²⁺-release events.

Effects of ß-adrenergic stimulation
ß-Adrenergic stimulation of cardiac myocytes has been found to result in a significant increase in both LCC and RyR Ca²⁺ influx, through a process involving channel phosphorylation mediated by the cAMP/PKA pathway (20,22,54). Under the effects of ß-adrenergic stimulation, active single-channel sweeps have shown that the increase in Ca²⁺ influx arises from a combined effect of elevated levels of open probability and availability, and is less likely to be a result of elevated unitary current values (23,43).

In this part of our study, we analyze the role of ß-adrenergic stimulation in Ca²⁺ signaling by introducing reaction schemes and parameters (Fig. 2, A and D, and Tables 1 and 2) relevant to channels under phosphorylated conditions. We simulate cases where the cleft molecular stoichiometry remains the same (five LCCs and 25 RyRs) while individual channels are either phosphorylated or nonphosphorylated, as specified. For example, suppose that persistent ß-adrenergic stimulation of the myocyte gives rise to phosphorylation of 20% of the channels. This implies that one LCC and five RyRs in the cleft will be represented by reaction schemes relevant to the phosphorylated case, whereas the other four LCCs and 20 RyRs will be governed by models for nonphosphorylated channels. Our ultimate goal is to describe the impact of channel phosphorylation on Ca²⁺ signaling. Consequently, we do not distinguish between the various ß-adrenergic receptor subtypes.

Gating of phosphorylated RyRs is a complex and controversial process that is yet to be fully understood. In the absence of consensus about the exact dynamics of phosphorylated RyRs, we consider a simple model to describe the activity of phosphorylated RyRs. Rather than making modifications to the model in Saftenku et al., which we use to describe nonphosphorylated RyRs, we choose to introduce a new reaction scheme (Fig. 2 D) similar in principle to that of the LCC. The phosphorylated RyR gates from a closed state (C.p) into an open state (O.p), whereas Ca²⁺ generation is mediated by transitions of the open state back to itself. We adopt reaction parameters (Table 2) from values reported in a study by Uehara et al. on the ligand-induced gating kinetics of phosphorylated RyR channels from canine myocytes (23). To describe the termination of SR Ca²⁺ release, we append a Ca²⁺-dependent inactivated state to the model (Fig. 2 D). Since we do not know how phosphorylation may alter the Ca²⁺ sensitivity of RyRs, we begin the phosphorylated RyR in an already-activated state. Studies by Reiken et al. (22) have correlated progressive PKA phosphorylation of RyR to cardiac dysfunction. We want this model to capture interesting effects of phosphorylation that may present a potential indicator for heart disease. We bear in mind that we have made major approximations in this part of our study, such that this set of results may broadly describe effects of phosphorylation but will not reveal other intricate Ca²⁺-release properties in all accuracy.

View larger version (38K): [in this window] [in a new window]   FIGURE 3  (A) Single-trial data on LCC flux, gating, and inactivation for the entire simulation duration (250 ms, left panels) and in higher time resolution (first ms, right panels). Data are obtained from a cleft model of Dcleft = 12 nm and Lcleft = 100 nm, containing five LCCs and 25 RyRs. No RyR activation or LCC voltage inactivation is observed in the first ms. The top traces track LCC flux, whereas the remaining traces track LCC activity in specific reaction states per time point. The top right trace also depicts the precise number of Ca2+ ions present in the cleft for the first ms, and when compared to the total LCC flux (of >10 times greater magnitude), can demonstrate that diffusion-mediated Ca2+ efflux occurs rapidly. (B) Representative visualization of activity (0.1 ms per frame) in the model described in Fig. 3 A. The TT membrane (top, shaded) bears five LCCs (green) that have all been Ca2+- or voltage-inactivated at 14 ms into the simulation. The SR membrane (bottom, unshaded) bears 25 RyRs in the closed (yellow) or open (red) states. Free Ca2+ ions are depicted as red spheres. This frame sequence illustrates a propagating pattern of RyR activation that is evident by following the open RyRs as time progresses.

Effects of Ca²⁺ buffering
The association and dissociation of free Ca²⁺ ions with endogenous Ca²⁺ buffers can be one of the determinants of Ca²⁺ spark properties. The majority of these buffers are immobile intracellular molecules with the potential to immobilize Ca²⁺ and momentarily dampen the effects of localized Ca²⁺ elevation. It has been estimated that >90% of Ca²⁺ ions entering a subcellular microdomain are able to bind rapidly to a range of immobile endogenous buffers constituting a variety of Ca²⁺-binding properties. For example, calmodulin and troponin C possess specific Ca²⁺-binding sites, whereas the TT and SR phospholipid membranes inherently possess high Ca²⁺-binding capacity but with low affinity (55).

To examine the effects of buffering on the local Ca²⁺ response, we model a few scenarios of buffering by incorporating immobile, membrane-bound Ca²⁺-binding sites on the top and bottom surfaces (TT and SR membranes). We compare the consequences of high and low densities of buffers and of buffers with high and low binding affinities, based on previously studied parameters (55).

The MCell simulation environment
MCell uses Monte Carlo algorithms (by pseudorandom number generation) to stochastically describe 3D Brownian random-walk diffusion and chemical reaction kinetics in complex spatial environments (www.mcell.cnl.salk.edu, (56)). Reaction transition probabilities are calculated according to user-specified rate constants, and compared to the value of a generated random number to decide the succeeding reaction state. The user defines the model using a special model description language, which MCell will parse to create corresponding C++ objects and simulations according to user instructions. To model a system using MCell, it is necessary to specify 1) geometry of the subcellular ultrastucture; 2) diffusion constants of diffusing ligands; 3) positions of effector sites that interact with ligands; 4) chemical reaction mechanisms and kinetic rate constants governing the system; and 5) an appropriate time step and the number of Monte Carlo time steps or iterations to simulate.

Simulations were performed on a dual Xeon 1.0 GHz workstation running the Hummingbird Exceed 8.0 X-server. It took 3 min of computer time to simulate 100 ms of real time. To speed up simulations, we used the MCell runtime optimization method of 3D spatial partitioning. Spatial partitions are transparent planes that the user places to create subvolumes in the modeled space, thus reducing the computation power required to track the movements of individual molecules to each subvolume. In this manner, computing speed is less dependent on microdomain complexity. MCell allows the user to export simulation results into visualization data formats for a variety of graphic tools. In this study, 3D images were rendered with IBM Data Explorer (www.opendx.org) using the companion visualization package DReAMM (www.mcell.psc.edu/DReAMM).

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES

 
Having established a Monte Carlo model of a dyadic cleft, one can try to investigate the properties of Ca²⁺ release within the cleft following a depolarization protocol representative of the voltage changes in an AP. Deterministic and partially deterministic models of myocyte Ca²⁺ release (4,7,57) have tracked LCC mediated responses to changes in membrane potential during the course of an AP using, e.g., Hodgkin-Huxley-type equations. Experimentally and theoretically, it has been determined that at the initiation of an AP, the cardiac ventricular sarcolemmal membrane experiences a sharp depolarization from around –80 mV to +20 mV followed by a gradual repolarization. Moreover, for more than three-quarters of an AP, the cardiac ventricular sarcolemmal membrane resides in the near-zero to positive potential range, implying that membrane potential does not change drastically much of this time, and that approximating the voltage dependence of LCC parameters will not critically affect dyadic cleft [Ca²⁺]. This assumption is further supported by findings that Ca²⁺-mediated inactivation is a more significant mode of LCC inactivation than voltage inactivation (58). In addition, although LCC currents play a major role in activating Ca²⁺ release from the jSR, they do not contribute significantly to the subspace [Ca²⁺] once SR release has been triggered. Thus Ca²⁺ spark duration and magnitude is relatively independent of the triggering sparklet (37). These considerations have prompted us to suggest that a simplified LCC simulation reflecting a voltage-clamp scenario can be a sufficiently accurate representation of trigger Ca²⁺ influx. An added benefit of this assumption is the possibility of defining the voltage-dependent parameters used in the model (unitary current, the rate of voltage-inactivation, and channel gating kinetics) to be equivalent to those experimentally determined under actual voltage-clamp conditions (0 mV, +10 mV, and +20 mV (37,44,43)). We thus assume throughout this study that the TT membrane is voltage-clamped at the beginning of each simulation to +10 mV and remains so for the entire response duration.

To confirm that LCC inactivation is primarily Ca²⁺-mediated, we tracked the LCC behavior in a dyadic cleft of L_cleft = 100 nm and D_cleft = 12 nm, containing five LCCs and 25 RyRs. At any time point, each of the five LCCs may exist in the open (conducting), closed, Ca²⁺-inactivated, or voltage-inactivated states. Fig. 3 A shows representative single-trial data on LCC flux, gating, and inactivation for the entire simulation duration (left panels), and, in greater detail, events that take place in the first millisecond (right panels). Note that in the initial phase of membrane excitation none of the RyRs have yet been activated by the trigger Ca²⁺. Since Ca²⁺ and voltage inactivation are known to be independent processes (46), here we assume that Ca²⁺-inactivated LCCs can also be voltage-inactivated. These traces show that Ca²⁺ inactivation of some of the LCCs can occur as early as <1 ms after the onset of the AP, whereas maximum Ca²⁺ inactivation (involving all five LCCs in the cleft) is reached within the first 20 ms. Voltage inactivation, on the other hand, gradually sets in later at 40 ms after the onset of the AP. Overall, the LCC opening events are distributed stochastically, with some instances of "bursts" representing a few openings in rapid succession. Note that although up to 650 Ca²⁺ ions are generated into the cleft in the first millisecond, the corresponding maximum Ca²⁺ level in the cleft amounts to only 50 ions. It is apparent that a vast majority of Ca²⁺ ions released into the modeled cleft are expected to diffuse rapidly into the surrounding cytosolic space. Fig. 3 B is a visual representation of sequential events in the cleft during the 14th ms, at a resolution of 0.1 ms per frame. At this time, all five LCCs (smaller hexagons, top surface) have been either Ca²⁺- or voltage-inactivated. The influx of Ca²⁺ ions (depicted as red spheres) in the cleft is observed to correspond to the spatial propagation of RyR (larger hexagons, bottom surface) activity. A more detailed view of Ca²⁺ influx over the first 5 ms is given in the Supplementary Material.

To determine the Ca²⁺ response within a single cleft after the onset of depolarization, we determined the number of Ca²⁺ ions within a standard cleft (L_cleft = 100 nm, D_cleft = 12 nm) as a function of time (Fig. 4 A). The single-cleft Ca²⁺ response (inset) shows that the number of Ca²⁺ ions can vary from zero to 60 ions over the course of a simulation, whereas the peak of the Ca²⁺ response averaged from 100 independent trials is only 23 ions. Thus, the average, predicted response is clearly far from deterministic and is subject to substantial stochastic variations. For the case of the standard cleft, the predicted Ca²⁺ response reaches a peak after 9.75 ms, gradually declining to basal levels by the end of 200 ms. This model reproduces critical features of experimentally observed Ca²⁺ sparks. In particular, the time to peak for the Ca²⁺ response is very similar to that previously recorded for Ca²⁺ sparks in mouse ventricular myocytes (10 ms) using confocal image analysis (59). At each given time, low numbers of ions in the cleft will translate into significant magnitudes of [Ca²⁺] due to the small cleft volume (1.2 x 10^–13 µL). The peak of the average response, for instance, is 300 µM. This value is significantly higher than the estimates presented in some previous deterministic and semistochastic models (4), but is consistent with other estimates (3). It is instructive to consider that a single Ca²⁺ ion within the assumed dyadic cleft volume would be equivalent to 13 µM, which is also above the range given in many modeling and experimental estimates. The significance of predicting these relatively high levels of Ca²⁺ in dyadic spaces will be discussed below.

View larger version (36K): [in this window] [in a new window]   FIGURE 4  The Ca2+ response after the onset of depolarization as a function of time. Each simulated cleft consists of five LCCs and 25 RyRs, and a standard Lcleft of 100 nm. (A) We first observe the response in a cleft with the assumed normal membrane separation of Dcleft = 12 nm, averaged from 100 independent trials. The averaged number of Ca2+ ions in the cleft (scatter plot) is fit to an exponential function (smooth curve) consisting of one rise-time constant and 1 decay-time constant. The predicted rise time of 9.75 ms is very close to that previously measured in Ca2+ sparks using rat ventricular myocytes. The inset shows results from two separate, single trials depicting the stochasticity of Ca2+ release. (B) Dcleft is progressively varied from 9–40 nm to study the effects of geometry changes in the cleft. After the onset of each AP, the Ca2+ level reaches a peak after 5–50 ms and gradually declines after 150 ms. The model demonstrates how, e.g., doubling of Dcleft from 12 nm to 24 nm, will result in a markedly reduced [Ca2+]. Further increasing the intermembrane distance to 40 nm produces a fluctuating trace of extremely low [Ca2+], a condition which might lead to altered contractile properties of the myocyte. The inset shows the dependence of rise time and peak [Ca2+] on Dcleft (x axis not to scale). This model corresponds qualitatively to results from earlier deterministic models of Ca2+ release.

Since the size of the dyadic cleft can also vary in terms of the dimensions of the TT and SR membranes (in addition to the gap separating them), we have also analyzed the dependence of Ca²⁺ response on the larger, x,y dimensions of the cleft (L_cleft). In particular, square membrane portions of 100–500-nm dimensions were analyzed for a constant D_cleft of 12 nm. Interestingly, the amplitude of Ca²⁺ response showed very little sensitivity to L_cleft in the 300–500 nm range, whereas the response was increasingly more pronounced as L_cleft was varied from 200 to 100 nm (Fig. 5 A). Overall, the amplitude of the response for the assumed D_cleft of 12 nm varied less than twofold. The rise-time variation was modest, changing marginally between 8.25 and 9.75 ms (Fig. 5 B). These results suggest that the dimensions of the TT and SR membranes might not be crucially important in determining the properties of a spark. Indeed, if one considers that changing L_cleft from 300 to 500 nm would result in an approximately threefold increase in the cleft volume and a 1.5-fold increase in the area of the "side walls" through which Ca²⁺ is assumed to escape, one sees that these geometry changes are comparatively similar to the increases in volume and escape area resulting from the increases in D_cleft from 9 to 25 nm and from 9 to 15 nm, respectively. However, the described change to L_cleft results in virtually no decrease in the peak [Ca²⁺], whereas the corresponding changes to D_cleft result in approximately six- and twofold decreases in the peak [Ca²⁺].

View larger version (22K): [in this window] [in a new window]   FIGURE 5  (A) [Ca2+] as a function of time in clefts containing five LCCs and 25 RyRs. Dcleft is fixed at 12 nm, whereas Lcleft ranges from 100–500 nm. (B) The dependence of peak [Ca2+] and rise time on Lcleft. The amplitudes of these parameters change modestly with Lcleft, in contrast to a higher dependence on Dcleft (Fig. 4 B, inset). These results suggest that membrane separation is likely to play a greater role in Ca2+ regulation than the lateral-cleft dimensions.

View larger version (21K): [in this window] [in a new window]   FIGURE 6  To investigate the relationship between the channel influx and the amount of Ca2+ present in the cleft, we obtain the integrals of these two parameters over time. (A) Correlations between the time integral of Ca2+ influx ((iLCC + iRyR)dt), and the sum total of Ca2+ residing in the cleft ( Ca2+dt), taken over the entire duration of each simulation. Lcleft ranges from 100 to 500 nm, whereas Dcleft is unchanged at 12 nm. (B) Dcleft ranges from 9 to 40 nm, whereas Lcleft is unchanged at 100 nm. Each line represents the correlation between channel influx and cleft Ca2+ for a particular geometry. When only Lcleft is changed, as in panel A, the magnitudes of the slopes exhibit a dependence on the cleft volumes. This is not observed when only Dcleft is varied, as in panel B. It is possible that greater membrane separation increases the cleft volume and its capacity to transiently store Ca2+, while compromising LCC/RyR coupling efficacy. The tradeoff between cleft volume and channel coupling can account for the resulting similar slopes for clefts of different volumes upon variation of Dcleft. (C) The time integral of Ca2+ influx through LCCs is compared with that through RyRs. No correlation is evident, suggesting that LCC currents serve only as a trigger regulating the onset of spark development, and that SR release is independent of LCC influx much of the time.

View larger version (28K): [in this window] [in a new window]   FIGURE 7  (A) Effects of varying geometry and channel density in models of Dcleft = 12 nm. Each CaRU is defined to be a cluster of five LCCs and 25 RyRs, and is denoted by a single asterisk in the figure. We first compare two models of similar dimensions (Lcleft = 200 nm; square x,y dimensions), with one containing four CaRUs (blue) and the other containing one CaRU (red). To investigate how local [Ca2+] will be affected when the spatial separation between clefts is removed, we also model a test case where Lcleft,1 = 100 nm and Lcleft,2 = 400 nm. This is equivalent to compacting multiple clefts adjacent to one another (see visualization in Fig. 1 C). (B) Comparison of models containing four CaRUs, across values of Dcleft ranging from 9 to 40 nm. The main figure depicts [Ca2+] as a function of time for the model of Lcleft = 200 nm, whereas the inset depicts the same for the model of Lcleft,1 = 100 nm and Lcleft,2 = 400 nm. Due to the rapid decay of Ca2+ response, the timescale for this figure is presented up to only 100 ms.

View larger version (16K): [in this window] [in a new window]   FIGURE 8  Effects of different levels of LCC and RyR phosphorylation in clefts containing five LCCs and 25 RyRs. The extent of channel phosphorylation was varied from 0% to 100% to investigate its influence on Ca2+ release (see Methods for detailed description). Phosphorylation of Ca2+ channels was found to progressively increase the amplitude of the peak response up to sixfold. One possible consequence is an alteration of myocyte contractile properties. Due to the rapid decay of Ca2+ response, the timescale for this figure is presented up to only 80 ms.

Establishing a fully stochastic model of Ca²⁺ regulation events in a dyadic cleft provides one with the advantage of being able to study the level of "noisiness" of the response. In particular, an important question relevant to the notion of local control of EC coupling is whether the localized Ca²⁺-release events are predicted to be too stochastic to exert any reasonable control of Ca²⁺ efflux from a dyadic space or of local contraction regulation. Another way this question can be put is "how distinct is the Ca²⁺ concentration in a single dyad from the average Ca²⁺ concentration produced by multiple sparks in multiple dyads that characterizes the whole cell level of contraction control?" If the deviation from the average Ca²⁺ level is strong, one can ask how many dyadic clefts should contribute to the local output to render the local response sufficiently close to the average to ensure robust contractility regulation. We have addressed this question by determining the number of Ca²⁺ ions in a standard dyadic cleft over 250 ms. It is evident (Fig. 4 A, inset) that single dyadic cleft responses can be significantly different in terms of both amplitude and duration of response from the average of multiple simulated dyadic-space responses (Fig. 4 A). It follows that to achieve the level of control suggested by the average [Ca²⁺] response, the outputs from several dyadic spaces may need to be coupled to attain the final contraction output. To assess the noise reduction with the number of summed up dyadic spaces we measured the deviation of the responses averaged over a varying number of clefts from the average of response for 100 clefts according to the following formula:

(1)

where is a time step in a simulation. We found that for the normal assumed cleft geometry X_n is reduced on average by 30%, if responses of 10 clefts are averaged, and by 50% for 20 clefts, compared to the average X_n value for a single cleft (Fig. 9). Thus the combined output can closely approach the true average spark response if 10–20 dyadic clefts contribute to it in a coordinated fashion. AP propagation can serve to synchronize the onset of the sparks so that the "averaging" is efficiently performed by the cell.

View larger version (20K): [in this window] [in a new window]   FIGURE 9  The stochasticity of Ca2+ signaling can be moderated if Ca2+ efflux from a relatively small number of clefts can be averaged in the downstream Ca2+-dependent events. The noise metric introduced in Eq. 1 is plotted as a function of the number of integrated dyadic clefts of different dimensions. Our estimates indicate that if the Ca2+ responses from 10–20 dyadic clefts are integrated, the noise in the Ca2+ output can be reduced by up to 50%.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION SUPPLEMENTARY MATERIAL ACKNOWLEDGEMENTS REFERENCES