Abstract
Many issues remain unresolved concerning how local, subcellular Ca^{2+} signals interact with bulk cellular concentrations to maintain homeostasis in health and disease. To aid in the interpretation of data obtained in quiescent ventricular myocytes, we present here a minimal whole cell model that accounts for both localized (subcellular) and global (cellular) aspects of Ca^{2+} signaling. Using a minimal formulation of the distribution of local [Ca^{2+}] associated with a large number of Ca^{2+}release sites, the model simulates both random spontaneous Ca^{2+} sparks and the changes in myoplasmic and sarcoplasmic reticulum (SR) [Ca^{2+}] that result from the balance between stochastic release and reuptake into the SR. Ca^{2+}release sites are composed of clusters of twostate ryanodine receptors (RyRs) that exhibit activation by local cytosolic [Ca^{2+}] but no inactivation or regulation by luminal Ca^{2+}. Decreasing RyR open probability in the model causes a decrease in aggregate release flux and an increase in SR [Ca^{2+}], regardless of whether RyR inhibition is mediated by a decrease in RyR open dwell time or an increase in RyR closed dwell time. The same balance of stochastic release and reuptake can be achieved, however, by either highfrequency/shortduration or lowfrequency/longduration Ca^{2+} sparks. The results are well correlated with recent experimental observations using pharmacological RyR inhibitors and clarify those aspects of the releasereuptake balance that are inherent to the coupling between local and global Ca^{2+} signals and those aspects that depend on molecularlevel details. The model of Ca^{2+} sparks and homeostasis presented here can be a useful tool for understanding changes in cardiac Ca^{2+ }release resulting from drugs, mutations, or acquired diseases.
 calcium signaling
 ryanodine receptor
 leak
 Markov chain
 tetracaine
intracellular Ca^{2+} is a ubiquitous biological signal that serves diverse functions in many cell types. In individual cells, information can be conveyed by both “global,” or cellwide, changes in [Ca^{2+}] and by “local,” subcellular Ca^{2+} signals. Local signals are frequently caused by the release of Ca^{2+} from intracellular stores, primarily the endoplasmic reticular (ER)/sarcoplasmic reticulum (SR). Local release events occur through closely packed clusters of release channels, inositol 1,4,5trisphosphate (IP_{3}) receptors or ryanodine receptors (RyRs), and are observable experimentally as Ca^{2+} sparks (7) or puffs (51). When one or several of the channels in a release site are open, the [Ca^{2+}] experienced by clustered channels is dramatically different from the [Ca^{2+}] in the bulk cytosol.
In cardiac myocytes and many other cell types, local and global Ca^{2+} signals are closely linked to one another. Local signals frequently form the building blocks from which global signals are built (32), and, conversely, changes in bulk [Ca^{2+}] in the myoplasm or SR can influence the frequency, amplitude, and kinetics of local events. In ventricular myocytes, for instance, propagating Ca^{2+} waves emerge when spontaneous Ca^{2+} sparks trigger additional sparks in a regenerative fashion (5). Similarly, changes in channel gating at the level of individual RyRs can immediately affect the production of local events and, over time, influence bulk myoplasmic and SR [Ca^{2+}]. An increase in RyR opening will cause a gradual decrease in SR [Ca^{2+}] (46), whereas inhibition of RyR opening, over time, will lead to elevated SR [Ca^{2+}] (21). This principle has been elegantly demonstrated by Lukyanenko and coworkers (35) in quiescent rat ventricular myocytes treated with caffeine (to sensitize RyRs) or tetracaine (to inhibit RyRs).
In heart cells, changes in Ca^{2+} signaling due to altered RyR activity are currently receiving considerable attention because of close links to disease (13, 48). In particular, catecholaminergic polymorphic ventricular tachycardia (CPVT), an inherited disorder associated with a dramatic increase in arrhythmia risk, results from mutations in either RyRs or calsequestrin, a SR Ca^{2+} buffer protein that associates with and modulates RyRs (20, 33). Experiments in vitro have shown that CPVTcausing mutations usually increase the open probability of the RyR, resulting in a hyperactive or “leaky” channel (12, 28, 31). Studies (1, 36) have also suggested that leaky RyRs are characteristic of several experimental heart failure models. Thus, a quantitative understanding of how changes in RyR gating influence local and global Ca^{2+} responses can provide insights into disease pathophysiology and can potentially suggest novel therapies.
The difference in spatial scales between local and global Ca^{2+} signals, however, creates significant challenges for the development of mechanistic mathematical models. In particular, gating of RyRs depends on both myoplasmic and SR [Ca^{2+}] (30), and concentrations within clusters during local events can be dramatically different from the bulk concentrations. In addition, because of the relatively small number of RyRs responsible for Ca^{2+} sparks (6), the stochastic gating of these channels must be considered when simulating local events. Previous studies (9, 27, 42, 44) have used Monte Carlo simulation methods to investigate the stochastic triggering of Ca^{2+} sparks, but these have generally treated myoplasmic and bulk SR [Ca^{2+}] as fixed boundary conditions. Conversely, modeling studies (4, 10, 40) focusing on cellular Ca^{2+} transients have usually used representations of SR Ca^{2+} release that do not account for the stochastic nature of the local events. Attempting to simulate Ca^{2+} signaling at both spatial scales simultaneously is a daunting prospect because the stochastic behavior of thousands of local events must be considered to determine the effects on the bulk concentrations. As a result, only a few studies (16, 17, 49, 50) have attempted to capture both phenomena.
In this report, we introduce a computationally efficient minimal model of the coupling between local and global Ca^{2+} signals in permeabilized ventricular myocytes. The model accounts for the random generation and termination of spontaneous Ca^{2+} sparks, the resulting changes in myoplasmic and SR [Ca^{2+}], and the feedback of these changes on spark frequency. We make few assumptions about the factors influencing RyR gating, which allows us to distinguish between those results that are inherent to the coupling between local and global signals and those that are specific to the RyR gating scheme. To validate the model, we considered experiments recently performed by Zima et al. (52), who showed that tetracaine, a RyR inhibitor, caused an initial suppression of Ca^{2+} sparks followed by an increase in SR [Ca^{2+}] and a partial recovery of spark frequency. Surprisingly, these authors found that prolonged exposure to tetracaine led to an increase in Ca^{2+} spark duration (see Fig. 1C in Ref. 52). The simulations presented here recapitulate these experimental results, suggesting that the observed increase in spark duration directly results from the interplay between RyR inhibition and the resulting changes in SR [Ca^{2+}]. More broadly, this model provides a powerful yet minimal framework for understanding how mutations, posttranslational modifications, or drugs can alter diastolic SR Ca^{2+} release in ventricular myocytes.
METHODS
Model formulation.
The minimal whole cell model of local and global Ca^{2+} responses developed here takes into account stochastic Ca^{2+}release site dynamics as well as the balance of release and reuptake fluxes leading to Ca^{2+} homeostasis in quiescent ventricular myocytes (Fig. 1). The model assumes that RyR Ca^{2+} channels are clustered on the ER/SR membrane in release sites composed of 10–30 channels. All channels in a given release site experience the same local [Ca^{2+}] (myoplasmic and SR), but these “domain” [Ca^{2+}] are heterogenous throughout the cell, i.e., different release sites experience different domain [Ca^{2+}]. Similar to previous work by Hinch and colleagues (16, 25, 26), we assume that when the number of open channels in a Ca^{2+}release site changes, the local [Ca^{2+}] rapidly equilibrates in a manner that balances the fluxes into and out of the spatially restricted domains. In our model formulation, a large number of stochastically gating Ca^{2+}release sites are coupled to the bulk myoplasmic and SR [Ca^{2+}] in a manner that allows spontaneous Ca^{2+} sparks to change the balance of Ca^{2+} release, or “leak,” and reuptake by sarco(endo)plasmic reticulum Ca^{2+}ATPase (SERCA) pumps. The bulk myoplasmic and SR [Ca^{2+}] determine the relationship between the number of open channels in a Ca^{2+}release site and the resulting domain [Ca^{2+}], and, consequently, changes in these bulk concentrations influence the dynamics of spontaneous sparks. This minimal yet realistic representation of bidirectional coupling between local and global aspects of Ca^{2+} handling is a novel aspect of our model formulation that has not been emphasized in previous work.
Ca^{2+}release site model.
In our model formulation, Ca^{2+}release sites are composed of N coupled Markov chains representing individual RyRs. For simplicity, we used a RyR model with two states, closed (C) and open (O), but the model formulation can be generalized for channel models with more states (18, 19). Each RyR opens at a rate that depends on the local myoplasmic domain (i.e., the diadic subspace) [Ca^{2+}], denoted by c_{myo}^{d}, and closes with a Ca^{2+}independent transition rate:
We assume that the local [Ca^{2+}] experienced by the Ca^{2+}regulatory site of each channel depends only on the number of open channels at the Ca^{2+} release site [N_{O}(t)] and the bulk Ca^{2+} concentrations (c_{myo} and c_{sr}), as described below. Because the channels are identical and indistinguishable, the state space for the N channel Ca^{2+}release site includes N + 1 states (0 ≤ N_{O} ≤ N), as follows:
Once the local [Ca^{2+}] (c_{myo}^{d,n}) that apply for any given number of open channels are specified, all of the statistical properties of the Ca^{2+}release site model can be determined using Eq. 3. In particular, the time evolution of the probability distribution of the number of open channels in the release site can be found by solving the following ordinary differential equation (ODE) initial value problem:
Myoplasmic and SR domain Ca^{2+}.
As shown schematically in Fig. 1, the domain [Ca^{2+}] for each release site (c_{myo}^{d} and c_{sr}^{d}) are coupled to the bulk compartments via the myoplasmic and SR fluxes (J_{myo} and J_{sr}) and coupled to one another through the release flux (J_{rel}). We assume these fluxes take the following form:
Recall that Eqs. 2 and 3 require the specification of the local [Ca^{2+}] that applies for any given number n of open channels (c_{myo}^{d,n}). Assuming the dynamics of domain Ca^{2+} are fast compared with the gating of RyRs, these domain [Ca^{2+}] are found by balancing the fluxes J_{rel}^{n} = J_{myo}^{n} and J_{sr}^{n} = J_{rel}^{n}. Solving these equations simultaneously for c_{myo}^{d,n} and c_{sr}^{d,n} yields the following:
The solid circles in Fig. 2 show the myoplasmic and SR domain concentrations given by Eqs. 8 and 9 as functions of the number of open channels (N_{O} = n). Note that as the number of open channels increases, c_{myo}^{d,n} increases and c_{sr}^{d,n} decreases, but always c_{myo}^{d,n} < c_{sr}^{d,n}. The open circles in Fig. 2 show how an increase in the bulk SR [Ca^{2+}] (c_{sr}) influences the relationship between the number of open channels in a Ca^{2+}release site and the resulting cytosolic (c_{myo}^{d,n}) and luminal (c_{sr}^{d,n}) domain [Ca^{2+}]. In particular, note that an increase in the bulk SR [Ca^{2+}] (c_{sr}) leads to an increase in cytosolic domain [Ca^{2+}] (c_{myo}^{d,n}) provided one or more RyRs are open (n ≥ 1 in Eq. 8).
Conventional Monte Carlo simulation.
As shown in Fig. 1, the bulk myoplasmic and SR [Ca^{2+}] are affected by the SERCA pump flux (J_{pump}), a passive leak from the SR to the myoplasm that is independent of release site activity (J_{leak}), as well as the fluxes from the myoplasmic domains to the bulk myoplasm (J_{myo}) and from the bulk SR to the SR domains (J_{sr}). In a conventional Monte Carlo simulation involving a large number of Ca^{2+}release sites, concentration balance equations for the myoplasmic and SR [Ca^{2+}] consistent with Fig. 1 are written as follows:
Accelerated simulation assuming a large population of release sites.
In our model formulation, we avoided Monte Carlo simulation of a large number of Markov chains by assuming that the number of release sites (and associated domains) is large enough that the probability distribution of release site states is well approximated by π(t) solving Eq. 4. In a naive application of this approach, Eqs. 11 and 12 would be integrated simultaneously with Eq. 4 with the substitution of π_{n} for f_{n} in Eqs. 13 and 14. Unfortunately, this approach is invalid because the assumption of rapid Ca^{2+} domain formation and collapse is a singular limit of the corresponding whole cell model formulation in which ODEs are used to solve for the dynamics of domain Ca^{2+} (see Supplemental Material).^{1} Instead, our model formulation is based on the following numerical solution of concentration balance equations for the total myoplasmic (ĉ_{myo}) and SR (ĉ_{sr}) [Ca^{2+}]:
Because the experimental observations that are of primary relevance to this report involve permeabilized cells (52), the plasma membrane flux (J_{pm}) was chosen to be as follows:
Summary and significance of the model.
Although minimal in nature, the whole cell model of local and global Ca^{2+} signaling that is the focus of this report accounts for the changes in myoplasmic and SR [Ca^{2+}] mediated by the balance of stochastic release and reuptake by the SR and the feedback of myoplasmic and SR [Ca^{2+}] on spark frequency. As discussed in the Introduction, previously published models of Ca^{2+} signaling in cardiac myocytes that include the stochastic release of SR Ca^{2+} either have not included bidirectional coupling between local Ca^{2+} release and global Ca^{2+} homeostasis or, because of the computational challenge of the required Monte Carlo simulations, have not emphasized the phenomenon. To our knowledge, this is the first systematic modeling study of the relationship between RyR kinetics, spontaneous and stochastic release of SR Ca^{2+}, and the resulting balance of bulk [Ca^{2+}] in permeabilized ventricular myocytes. It is also the first model of Ca^{2+} sparks and homeostasis that bypasses Monte Carlo simulation by assuming both a large number of Ca^{2+}release sites and rapid Ca^{2+} domain dynamics, resulting in a minimal formulation that facilitates parameter studies.
The minimal whole cell model of local and global Ca^{2+} signaling that is the focus of this report includes N + 3 ODEs and several algebraic relations. Two ODEs are concentration balance equations for the total myoplasmic (ĉ_{myo}) and SR (ĉ_{sr}) [Ca^{2+}] (Eqs. 15 and 16). The additional N + 1 ODEs (Eq. 4) account for the dynamics of a large number of Ca^{2+}release sites, each composed of N twostate RyRs (Eqs. 1 and 2). Algebraic relations include the fluxes (Eq. 17 and Eqs. 22–24) that appear in the concentration balance equations as well as the assumed relationship between myoplasmic (c_{myo}^{d,n}) and SR (c_{sr}^{d,n}) domain [Ca^{2+}] and the number of open channels (Eqs. 8–10). Note that the fluxes J_{rel}^{T}, J_{pump}, and J_{leak} are functions of c_{myo} and c_{sr}, which are functions of ĉ_{myo}, ĉ_{sr}, and π. The algebraic relationship between these quantities is found by inverting Eqs. 18 and 19 after substitution of Eqs. 8, 9, 20, and 21 (see Supplemental Material).
Although our model formulation assumes a large population of Ca^{2+}release sites, we do not have to specify a precise number. To see this, note that the domain concentrations c_{myo}^{d,n} and c_{sr}^{d,n} do not depend on the number of release sites in the cell (M) when the rate constants are defined by v_{rel} = v_{rel}^{T}/M, v_{myo} = v_{myo}^{T}/M, and v_{sr} = v_{sr}^{T}/M (Eqs. 8 and 9). Because the rate v_{rel}^{T} that appears in Eq. 17 does not correspond to release through one release site but rather the entire population, it is convenient to specify c_{myo}^{d,n} and c_{sr}^{d,n} using Eqs. 8–10 with the replacement of v_{rel}^{T}, v_{myo}^{T}, and v_{sr}^{T} for v_{rel}, v_{myo}, and v_{sr}.
The minimal model of local and global Ca^{2+} signaling includes 14 parameters, far fewer than most mathematical models of Ca^{2+} handling in cardiac myocytes (see the Supplemental Table in the Supplemental Material). Some parameters [such as the effective volume ratios λ_{sr}, Λ_{sr}^{d}, and Λ_{myo}^{d} and the SERCA pump maximum rate (v_{pump}) and dissociation constant (k_{pump})] were either chosen to be consistent with previous work (49, 50) or do not require extensive consideration because model responses to changes in these parameters are obvious and intuitive. Because the ventricular myocyte is assumed to be permeabilized, the precise value of the parameter k_{pm} is unimportant so long as there is rapid equilibration of bulk myoplasmic Ca^{2+} with the extracellular [Ca^{2+}] (c_{ext}). The assumed number of RyRs in each release site (N = 10) was chosen to be consistent with estimates of the number of channels activated during a Ca^{2+} spark (7, 15, 43). This is a smaller number of RyRs than previously reported in electron microscopic studies performed a decade ago (11, 34) but is consistent with more recent estimates based on superresolution optical techniques and threedimensional electron tomography (2, 23). The most important of the model parameters [the kinetic parameters for the stochastic gating of the twostate RyR (k_{co} and k_{oc}) and the rate constants for Ca^{2+} release v_{rel}^{T}, myoplasmic domain collapse v_{myo}^{T}, and luminal domain recovery v_{sr}^{T}] are more difficult to constrain, and, consequently, these parameters are the focus of numerous sensitivity studies (see below).
The following aspects of the model behavior suggest that our standard parameter set is physiologically realistic. At 100 nM cytosolic [Ca^{2+}] (c_{myo}), the average duration of a spontaneous Ca^{2+}release event is on the order of 20 ms, similar to the observed rise time of Ca^{2+} sparks (6, 7). The Ca^{2+} spark rate with c_{myo} = 100 nM is 0.043 sparks/s per release site (∼1 spark every 23 s). Assuming 20,000 release sites in a ventricular myocyte, this corresponds to 860 sparks/s per cell, that is, 86 sparks/s in a fast twodimensional confocal frame scan that samples 10% of the cell volume. This value is consistent with an experimental study (3) performed in intact cells that reported spontaneous spark rates of 1–4 × 10^{−5} μm^{−2}·ms^{−1}, which corresponds to 30–120 sparks/s assuming a crosssectional area of 100 × 30 μm. Consistent with experiment, an increase in myoplasmic [Ca^{2+}] in the model leads to an increase in the spontaneous Ca^{2+} spark rate.
RESULTS
RyR open probability and spontaneous Ca^{2+} sparks.
The minimal model of local and global Ca^{2+} signaling that is the focus of this report simulates stochastic Ca^{2+} release by clusters of RyRs and the resulting whole cell Ca^{2+} homeostasis in quiescent ventricular myocytes. The modeling formalism (described in Summary and significance of the model) was chosen to be as simple as possible while still provide mechanistic insights into the perturbation of SR Ca^{2+} leak that results from pharmacological agents, mutations, or posttranslational modifications of the RyR, as may occur in disease states. For example, tetracaine, a potent local anesthetic that allosterically blocks Ca^{2+}release channels, reduces the open probability of RyRs in planar lipid bilayer experiments (52) by increasing the mean closed dwell time of channels (21). Because the mean closed time of the twostate RyR model is given by τ_{C} = 1/k_{co}(c_{myo}^{d})^{2} (Eq. 1), we simulated the application of tetracaine to permeabilized ventricular myocytes by decreasing the rate constant k_{co}, which influences the stochastic dynamics of the Ca^{2+}release sites (Eq. 2) in the minimal whole cell model. We wanted to understand how the simulated application of tetracaine influences the dynamics of Ca^{2+} sparks and homeostasis in the permeabilized ventricular myocyte model.
Figure 3 shows a summary of 60 numerical calculations of the stationary dynamics of the minimal whole cell model performed using different values of the RyR Ca^{2+} activation rate constant k_{co}. The circles show the result of two particular simulations: one corresponding to the standard parameter values (k_{co} = 4.5 μM^{−2}·s^{−1}; solid circles) and the other corresponding to the simulated application of tetracaine (k_{co} = 0.5 μM^{−2}·s^{−1}; open circles). In the latter case, the value of k_{co} was chosen so that the single channel probablity (P_{open}), given by the following:
Figure 3A shows that the simulated application of tetracaine led to increased bulk SR [Ca^{2+}] in the whole cell model (c_{sr} = 342 to 1,112 μM; compare solid and open circles). As expected, steadystate bulk SR [Ca^{2+}] increased as RyR open probability decreased, due to a decrease in the total release flux (Eq. 17). Because SERCA pump flux (Eq. 23) is independent of bulk SR [Ca^{2+}], maximum bulk SR [Ca^{2+}] (c_{sr}) asymptotically approaches 2.5 mM when association rate constant k_{co} is very small. When a nonspecific passive leak was not included in the model, c_{sr} increased further (not shown). Results similar to those shown in Fig. 3A can be obtained without a passive leak by extending the SERCA pump model to include both forward and reverse modes (39).
Figure 3B shows that during the simulated application of tetracaine, the fraction of open channels in a randomly sampled release site (f_{O}) was reduced by 79%, less than the reduction in the single channel open probability given by Eq. 25 (88%). This result indicates that elevated bulk SR [Ca^{2+}] and the interaction between RyRs combine to attenuate the decrease in channel activity occurring during the simulated application of tetracaine. That is, increased SR [Ca^{2+}] increases the driving force during stochastic Ca^{2+}release events and elevates myoplasmic domain [Ca^{2+}] (compare solid and open circles in Fig. 2B). The fraction of open channels can be calculated as f_{O} = E[N_{O}]/N, where E[N_{O}] is the expected value of the number of open channels in a randomly sampled release site:
With spark initiation defined as a release site reaching a threshold number of open channels (N_{O} = 4 to 5 transition) and spark termination defined as all channels closing (N_{O} = 1 to 0 transition), spark frequency and mean duration were calculated using the matrix analytic method described in Ref. 18. The solid and open circles in Fig. 3, C and D, show that the simulated application of tetracaine decreased Ca^{2+} spark frequency but increased mean Ca^{2+} spark duration, consistent with experimental observations (52). While bulk SR [Ca^{2+}] (c_{sr}) and spark frequency are monotone functions of RyR open probability (decreasing and increasing, respectively), mean spark duration is a biphasic function of RyR open probability (first increasing and then decreasing).
Spark frequency and duration upon the application of tetracaine.
Figure 4A, left, shows representative Ca^{2+}release events exhibited by a Ca^{2+}release site in the standard whole cell simulation. Spontaneous Ca^{2+} sparks were simulated using Gillespie's method (14). The parameters used correspond to the solid circles shown in Fig. 3 (k_{co} = 4.5 μM^{−2}·s^{−1}) and resulted in robust Ca^{2+} sparks (score = 0.51) for the bulk concentrations (c_{myo} and c_{sr}) of the equilibrated whole cell model. Note the high frequency of spontaneous release events, including five sparks (N_{O} ≥ 5) and a large number of smaller release events. These small release events, termed “Ca^{2+ }quarks” (32), would not be detectable with standard confocal microscopy and would therefore contribute to “invisible” SR Ca^{2+} leak (43). Figure 4A, right, shows an expanded version of the first spark (asterisk in Fig. 4A, left), which has a duration (17.5 ms) close to the mean spark duration with these standard parameters (18.3 ms; solid circle, Fig. 3D).
Figure 4B shows representative Ca^{2+}release events during the simulated addition of tetracaine. The parameters used correspond to the open circles shown in Fig. 3 (k_{co} = 0.5 μM^{−2}·s^{−1}). While robust Ca^{2+} sparks were observed (score = 0.51), the simulated application of tetracaine significantly reduced spark frequency; one spark and three quarks were observed. While the mean spark duration was 26.8 ms (open circle, Fig. 3D), Fig. 4B, right, shows an expanded version of the observed spark, which was over 90 ms in duration. Consistent with experimental observations (52), such longduration sparks are not infrequent during the simulated addition of tetracaine, despite the fact that they almost never occurred with the standard parameter set (see below).
To confirm that this decreased spark frequency and increased mean spark duration are due to overloading of bulk SR [Ca^{2+}], Fig. 4C shows a control simulation using the single channel RyR parameters shown in Fig. 4B (k_{co} = 0.5 μM^{−2}·s^{−1}) with bulk SR [Ca^{2+}] “clamped” at the value shown in Fig. 4A (c_{sr} = 342 μM). The resulting simulation showed only a few release events, none with more than three channels open (score = 0.11). We conclude that the overloading of SR [Ca^{2+}] that occurs when RyR open probability is decreased is required for the presence of prolonged sparks in the whole cell model.
As mentioned above, the simulated addition of tetracaine resulted in Ca^{2+} sparks whose duration tended to be longer than that observed with the standard parameters (compare Fig. 4, A, right, and B, right). To further quantify this effect, Fig. 5A shows the numerically calculated distribution of spark durations for the standard (solid line) and tetracaine (dashed line) parameter sets. While the mode of these distributions was nearly identical (standard: 4.8 ms and tetracaine: 4.6 ms), in the case of tetracaine the distribution extended further to the right, consistent with the higher probability of long sparks (Fig. 4B, right). Integrating the results shown in Fig. 5A led to cumulative probability distributions (Fig. 5B) that showed that 21.5% of the sparks in the tetracaine simulations but only 9.4% of the sparks in the standard simulations were longer than 40 ms (compare solid and open circles).
Magnitude of Ca^{2+} release due to spontaneous sparks.
As discussed above, the SR [Ca^{2+}] overload induced by tetracaine led to higher myoplasmic domain [Ca^{2+}] (compare open and solid circles in Fig. 2A), higher SR domain [Ca^{2+}] (Fig. 2B), and higher release flux (J_{rel}^{T}, Eq. 17) for any given number of open channels. The resulting changes in the dynamics of Ca^{2+}mediated coupling of RyRs led to a decrease in spark frequency and an increase in spark duration (Fig. 3, C and D). Nevertheless, the solid and open circles in Fig. 5C show that the application of tetracaine decreases the aggregate release flux in the whole cell model. In fact, the solid line in Fig. 5C shows that the aggregate release flux is a monotone increasing function of RyR open probability, despite the fact that the SR [Ca^{2+}] is monotone decreasing (Fig. 3A).
During experimental observations of spontaneous Ca^{2+} release, smallamplitude events may not be detectable. Thus, it is of interest to dissect the aggregate release flux of the whole cell model to determine the fraction of spontaneous release that occurs via release sites that have few open channels. The dotted and dashed lines in Fig. 5C show that the release flux mediated by release sites with one or two open channels is a monotone increasing function of the RyR Ca^{2+} activation rate constant k_{co}. Figure 5D shows aggregate flux J_{rel}^{T} jointly distributed with the number of open channels for both the tetracaine and standard parameter sets. Both distributions were bimodal; a peak was observed at N_{O} = 1 as well as N_{O} = 6 or 7. Tetracaine suppressed the proportion of Ca^{2+} released through sites with seven or fewer open channels (N_{O} ≤ 7), whereas release mediated by sites with eight or more open channels (N_{O} ≥ 8) increased slightly in the presence of tetracaine. These observations are consistent with the increased probability of longduration Ca^{2+} sparks observed upon the application of tetracaine (Fig. 5A).
Because the detectability of sparks recorded with fluorescent dyes is primarily determined by spark amplitude (i.e., integrated Ca^{2+} release) (8, 41), Fig. 6A shows a summary of Monte Carlo simulations analyzing how spark amplitude and duration are jointly influenced by the simulated application of tetracaine (compare solid and open circles). Here, spark events were defined as beginning with a N_{O} = 0 → 1 transition, whereas spark amplitude was the integrated stochastic release flux (Eq. 5) before spark termination via a N_{O} = 1→ 0 transition. Figure 5, B and C, shows the cumulative probability distributions of spark amplitude and duration, respectively. Sparks had both larger amplitude and extended duration when tetracaine was applied (consistent with Fig. 5A). Nevertheless, the decrease in spark frequency in the presence of tetracaine led to an overall decrease in the aggregate release flux discussed above (Fig. 5C).
Figure 6D shows the percentage of undetected spark events (and the “hidden” Ca^{2+}release flux mediated by undetected sparks) as a function of a detection threshold on spark amplitude. The vertical dashed lines in Fig. 6, B and D, indicate a sensitive detection threshold equivalent to the amount of Ca^{2+} released by an average quark; in the standard simulation, this is a single channel release event with a duration of 1.1 ms. With this sensitive detection threshold, 36% of release events are not observed; most of these release events are quarks, brief single channel openings through which <1% of the stochastic Ca^{2+} release occurs. Because the simulated application of tetracaine led to increased SR load and greater release flux for any given number of open RyRs (Fig. 2), the tetracaine condition led to fewer hidden events (15%) and decreased hidden release (<0.1%). The solid and dashed lines in Fig. 6D show that the percentage of hidden release events and hidden release flux were both increasing functions of detection threshold. For the range of possible detection thresholds shown, the percentage of hidden events decreased by two to threefold upon the application of tetracaine.
Transient effects upon the application of tetracaine.
While the above simulations focused on steadystate dynamics of the whole cell model, Fig. 7 shows transient effects upon bulk SR [Ca^{2+}] (c_{sr}), mean spark duration, and spark frequency that occurred during the simulated application and washout of tetracaine. Consistent with experimental observations (21, 52), the initial application of tetracaine caused spark frequency to decrease; the mean spark duration during this phase (5.0 ms) was much shorter than the baseline value (18.3 ms). However, this reduced spontaneous Ca^{2+} release caused a slow increase in SR load that ultimately increased the mean spark duration to 26.8 ms, consistent with the steadystate results (Fig. 3D).
Figure 7 also shows that upon the simulated washout of tetracaine (right arrow), there was a transient increase in spark frequency (maximum of ∼4 times the baseline value) and a rapid depletion of SR [Ca^{2+}] from elevated to baseline values. For a short period of time, the mean spark duration was quite large (see asterisk); however, the value attained is not relevant because it is greater than duration of the phase itself (400 ms). Shortly after this burst of spark activity, the mean spark duration returned to baseline.
Model parameters and Ca^{2+} homeostasis.
Figure 8 shows a summary of 2,500 calculations of the stationary dynamics of the whole cell model as a function of v_{sr}^{T} and v_{myo}^{T}. These parameters control the rate of “diffusion” or translocation of Ca^{2+} between the different cellular subspaces represented in the minimal model, for example, from the individual myoplasmic domains (diadic subspaces) to the cytoplasm and from the bulk SR to the luminal domains (network to junctional SR). These domain rate constants also influence the strength of the Ca^{2+}mediated coupling between RyRs and the extent of SR domain depletion during sparks (recall Eqs. 8–10). Because v_{sr}^{T} and v_{myo}^{T} are difficult to constrain via experiments, they are good choices for a parameter study designed to determine their effect on experimentally observable quantities such as steadystate bulk SR [Ca^{2+}] and mean spark frequency and duration.
Figure 8A shows that bulk SR [Ca^{2+}] is an increasing function of rate of myoplasmic domain collapse (v_{myo}^{T}) and a decreasing function of the rate of SR domain recovery (v_{sr}^{T}). When v_{myo}^{T} is very large or v_{sr}^{T} is very small, RyRs become decoupled such that most openings fail to trigger neighboring channels, and the reduced leak causes an increase in bulk SR [Ca^{2+}]. Figure 8, B and D, shows that the fraction of open channels and mean spark duration were much more sensitive to v_{myo}^{T} than to v_{sr}^{T}. When v_{myo}^{T} is small, sparks are extremely long because subspace [Ca^{2+}] remains elevated after RyRs close (cf. Refs. 27 and 37). Fast SR refilling alone is not sufficient to induce longduration Ca^{2+} sparks of the sort observed upon the simulated application of tetracaine (cf. Fig. 4B).
Robust Ca^{2+} sparks were observed in the whole cell model even when the domain rate constants were ranged over several orders of magnitude (Fig. 9A). This robust Ca^{2+} spark behavior is a consequence of the homeostatic changes in bulk SR [Ca^{2+}] accounted for in our model formulation; when these simulations were repeated with bulk myoplasmic and SR Ca^{2+} concentrations fixed at baseline values (c_{myo} = 0.1 μM and c_{sr} = 342 μM), the range of domain rate constants leading to Ca^{2+} sparks was considerably smaller (Fig. 9B).
Figure 10 shows stationary dynamics of the whole cell model as a function of bulk myoplasmic [Ca^{2+}] (c_{myo}) and maximum rate of the SERCA pump (v_{pump}), two parameters that can be easily manipulated in experiments. Bulk SR [Ca^{2+}] was (as expected) an increasing function of v_{pump}; however, the SR Ca^{2+} load was a biphasic function of c_{myo} (Fig. 10A). The fraction of open channels and mean spark duration were both increasing functions of SERCA pump activity (Fig. 10, B and D). Spark frequency was not sensitive to SERCA activity but was a rapidly increasing function of bulk myoplasmic [Ca^{2+}] (Fig. 10C). Spark duration was a biphasic function of c_{myo} (Fig. 10D), consistent with the biphasic effect of this parameter on SR load (Fig. 10A). Bulk myoplasmic [Ca^{2+}] (c_{myo}) significantly changed the functional dependence of the stationary dynamics of the whole cell model on the RyR Ca^{2+} activation rate constant (k_{co}). On the other hand, when the application of tetracaine was simulated by comparing k_{co} values that correspond to a decrease in RyR activity from f_{O} = 9.6 × 10^{−4} to 2.0 × 10^{−4} (Table 1), we observed increased SR load, decreased spark frequency, and increased spark duration for c_{myo} in the range of 0.1–0.5 μM.
RyR inhibition mechanism and spark duration.
We modeled the action of tetracaine as a decrease in the Ca^{2+} activation rate constant k_{co}, which reduces the open probability (Eq. 25) of the RyR model (Eq. 1) by increasing the mean closed dwell time, τ_{C} = 1/k_{co}(c_{myo}^{d})^{2}. However, the open probability of the RyR can also be reduced by increasing the rate constant k_{oc}, thereby decreasing the mean open dwell time (τ_{O} = 1/k_{oc}). Such a change would be analogous to the pharmacological action of the antiarrhythmic agent flecainide, which has been shown to reduce the dwell time of RyR open states (47).
Figure 11, C and D, shows Ca^{2+} spark frequency and mean spark duration as a function of the RyR kinetic constants k_{co} and k_{oc} when these spark statistics are well defined (score ≥ 0.2). The contours dividing the plane into the areas where sparks are present (gray) and absent (white) follow constant K = k_{oc}/k_{co}, indicating that sparks occur provided the single channel RyR open probability (Eq. 25) is neither too low or too high. If RyR parameters were changed from the standard values (asterisk), the resulting change in SR [Ca^{2+}] also depend only on RyR open probability (Fig. 11A). SR [Ca^{2+}] was too high for sparks in the top left white region and too low for sparks in the bottom right (cf. Fig. 3A). In contrast, spark frequency and mean spark duration strongly depended on whether a pharmacological perturbation of RyR kinetics is assumed to affect k_{co}, k_{oc}, or both (diamonds). A reduction in k_{co} (increase in τ_{C}), analogous to effect of tetracaine, led to fewer sparks but an increase in spark durations. Conversely, an increase in k_{oc} (decrease in τ_{O}), analogous to the effect of flecainide, caused more sparks with decreased mean duration. Examples of MonteCarlo simulations under these different conditions are shown in Fig. 12. Similarly, we can assume that lowdose caffeine leads to an increase k_{co} accompanied by a smaller decrease in k_{oc} (+ symbol) (29). The model predicted that this would cause a decrease in SR [Ca^{2+}], an increase in spark frequency, and little change in spark duration, roughly consistent with the results observed by Lukyanenko et al. (35).
These results illustrate a fundamental point about the interplay between local and global Ca^{2+} signals during pharmacological interventions designed to manipulate spontaneous cellular responses. Local (microscopic) aspects of cell response, such as spark frequency and duration, are highly dependent on the molecular details, that is, precisely how of the kinetics of RyR stochastic gating kinetics have been modified. Global (macroscopic) aspects of the cell response, such as steadystate bulk SR [Ca^{2+}], are less dependent on the kinetic details of the perturbation but remain determined by equilibrium quantities such as the dissociation constant for Ca^{2+} binding to the RyR.
DISCUSSION
This report presents a minimal whole cell model of local and global Ca^{2+} signals in quiescent ventricular myocytes. The modeling formalism accounts for the effect of random spontaneous Ca^{2+} sparks, changes in bulk myoplasmic and SR [Ca^{2+}] mediated by the balance of stochastic release and reuptake by the SR, and feedback of myoplasmic and SR [Ca^{2+}] on spark frequency. The functional organization of the model (Fig. 1) is similar to previously published Monte Carlo models of local control (17). The assumptions made here regarding the rapid equilibration of domain Ca^{2+} are similar to assumptions made in previously published local control models that represent the stochastic dynamics of a large number of Ca^{2+}release sites (16, 25), but these previous studies have not made a distinction between domain and bulk SR Ca^{2+} as done in our minimal whole cell model.
To our knowledge, this is the first theoretical study of the relationship between RyR kinetics, spontaneous and stochastic Ca^{2+} release, and the resulting balance of bulk myoplasmic and SR [Ca^{2+}] in quiescent ventricular myocytes. Because of the computational challenge of largescale simulations, a traditional Monte Carlo approach is not well suited to investigate these phenomena. The whole cell modeling approach introduced here bypasses Monte Carlo simulation by assuming a large number of Ca^{2+}release sites and rapid Ca^{2+} domain dynamics, resulting in a minimal formulation that facilitates parameter studies.
The minimal model of local and global Ca^{2+} responses in quiescent ventricular myocytes presented here is able to recapitulate recent experiments by Zima et al. (52) showing that tetracaine, an inhibitor of RyRs, causes a transient suppression of Ca^{2+} sparks that is followed by an increase in bulk SR [Ca^{2+}], partial recovery of spark frequency, and an increase in Ca^{2+} spark duration (Fig. 11). Conversely, if flecainide is assumed to decrease RyR mean open time while not affecting closed times, the model predicts that this drug will cause an increase in spark frequency and decrease in spark duration, similar to recent experimental observations (24). However, it should be noted that this study found no change in steadystate SR [Ca^{2+}], suggesting that the effects of flecainide on RyR gating are somewhat more complex.
More broadly, these examples illustrate how the model provides insights into the relationships between RyR gating, Ca^{2+} spark characteristics, and balance of bulk myoplasmic and SR [Ca^{2+}] (Figs. 3⇑–5). For example, simulations showing that mean spark duration is a biphasic function of the RyR Ca^{2+} activation rate constant k_{co} suggest that the increase in spark duration observed after the application of tetracaine may be concentration dependent (Fig. 3D). Despite the fact that spark duration is biphasic, a dissection of model responses suggests that the release flux is a monotone function of RyR open probability. The model also predicts that tetracaine suppresses the hidden flux mediated by Ca^{2+}release events below detection threshold more strongly than observable release events (Fig. 6D). Parameter experiments exploring the effect of different mechanisms of RyR inhibition (Table 1 and Fig. 12) indicate that the whole cell Ca^{2+} balance is largely determined by RyR open probability (a function of K = k_{oc}/k_{co}). Conversely, the frequency and duration of Ca^{2+} sparks are sensitive to RyR open and closed dwell times, which are determined by k_{co} and k_{oc} independently (see Fig. 11).
The Ca^{2+}release site model used here is quite simple and assumes “instantaneous mean field coupling” (38) of N = 10 twostate RyRs. When model simulations were performed with release sites composed of larger numbers of channels, qualitatively similar results were observed (supplemental Fig. S1). Our whole cell modeling approach can be applied with single channel models of arbitrary complexity, so long as the state space of the resulting Ca^{2+}release site model is not so large that integrating Eq. 4 is impractical. Our use of a twostate RyR in this first study of the relationship between Ca^{2+} sparks and homeostasis allows us to focus on aspects of the cellular response to perturbations that are likely to be fundamental and general. A more complicated RyR model would call into question the specific details of how we (arbitrarily) suppose pharmacological perturbations influence the kinetic constants that determine RyR stochastic gating, assumptions that could easily influence the whole cell response and our conclusions. We find it intriguing that this model of local and global Ca^{2+} responses in quiescent ventricular myocytes is able to reproduce changes in spark frequency and duration caused by the application of both tetracaine and flecainide, despite the fact that the RyR model used does not include regulatory processes such as Ca^{2+}dependent inactivation and/or sensitization by SR [Ca^{2+}]. The effect of these wellestablished aspects of RyR Ca^{2+} regulation is beyond the scope of this work.
The Ca^{2+} activation process in the RyR model is mediated by myoplasmic domain Ca^{2+} (c_{myo}^{d}, Eq. 9), and it would be possible to augment the model to include the sensitization of RyRs by SR domain Ca^{2+} (c_{sr}^{d}, Eq. 9). However, the assumption of instantaneous coupling of RyRs (i.e., rapid equilibration of myoplasmic and SR domain Ca^{2+}) may not work well when luminal Ca^{2+} plays a role in spark termination (27, 42). A computational study of the contribution of luminal RyR regulation to the bidirectional coupling of spontaneous Ca^{2+} release and Ca^{2+} homeostasis will likely require more subtle mathematical formulations, similar to the probability density and moment closure techniques that accelerate “local control” simulations of highgain graded Ca^{2+} release in voltageclamped cardiac myocytes (49, 50). These representations of heterogenous domain Ca^{2+} concentrations associated with a large number of Ca^{2+}release sites remain valid even when the dynamics of SR domain Ca^{2+} are not fast compared with channel gating. Indeed, the model presented here can be viewed as a reduction of the moment closure formulation (50) that is valid when SR Ca^{2+} domains rapidly equilibrate with myoplasmic domain and bulk SR Ca^{2+}, in which case there is a negligible variance of SR domain [Ca^{2+}] for any given Ca^{2+}release site state.
In preliminary work, we studied whole cell responses using RyR models that included both Ca^{2+} activation and inactivation by myoplasmic domain Ca^{2+}, for example, the following threestate model:
While the simulated spark duration histograms shown in Fig. 5A are unimodal, the experimentally observed spark duration histogram shows two peaks in the presence of tetracaine, suggesting two distinct populations of sparks (52). Zima and coworkers (52) suggested that the prolonged sparks occurred at release sites with highly interconnected junctional SR and high refilling rates. Our simulations shown in Fig. 8 address this idea by investigating the effects of a larger domain refilling rate (v_{sr}^{T}). Consistent with the hypothesis of Zima et al. (52), an increase in v_{sr}^{T} from 10 to 50 s^{−1} caused a slight increase in baseline Ca^{2+} spark duration as well as a greater percent increase upon the application of tetracaine (57% vs. 47%, not shown). This suggests that the long sparks observed upon the application of tetracaine (52) may indeed be associated with fast rather than slow SR refilling.
Finally, it is intriguing that robust Ca^{2+} sparks were observed in the whole cell model even when the domain rate constants were ranged over several orders of magnitude (Fig. 8). Because the range of domain rate constants leading to Ca^{2+} sparks is considerably smaller when bulk myoplasmic and SR [Ca^{2+}] are fixed at baseline values (compare Fig. 9, A and B), we conclude that this robust Ca^{2+} spark behavior is a consequence of homeostatic changes in bulk SR [Ca^{2+}] accounted for in our model formulation. The fact that the feedback of spontaneous SR leak on spark frequency extends the regime where spontaneous sparks occur is potentially physiologically relevant. Speaking teleologically, the homeostatic mechanisms appear to encourage SR leak mediated by Ca^{2+} sparks and discourage the alternatives: SR leak mediated by quarks or tonically active release sites. These observations underscore the importance of accounting for global Ca^{2+} balance in models of localized Ca^{2+}release events.
GRANTS
This work was supported in part by National Science Foundation (NSF) Grant 0443843 (to G. D. Smith) and a Howard Hughes Medical Institute grant through the Undergraduate Biological Sciences Education Program to The College of William and Mary. This work was performed in part using computational facilities at The College of William and Mary provided with the assistance of the NSF, Virginia Port Authority, Sun Microsystems, and Virginia's Commonwealth Technology Research Fund.
DISCLOSURES
No conflicts of interest, financial or otherwise, are declared by the author(s).
ACKNOWLEDGMENTS
Some of the results of this study have previously appeared in abstract form (22).
Footnotes

↵1 Supplemental Material for this article is available online at the American Journal of PhysiologyHeart and Circulatory Physiology website.
 Copyright © 2010 the American Physiological Society