Many issues remain unresolved concerning how local, subcellular Ca2+ 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 Ca2+ signaling. Using a minimal formulation of the distribution of local [Ca2+] associated with a large number of Ca2+-release sites, the model simulates both random spontaneous Ca2+ sparks and the changes in myoplasmic and sarcoplasmic reticulum (SR) [Ca2+] that result from the balance between stochastic release and reuptake into the SR. Ca2+-release sites are composed of clusters of two-state ryanodine receptors (RyRs) that exhibit activation by local cytosolic [Ca2+] but no inactivation or regulation by luminal Ca2+. Decreasing RyR open probability in the model causes a decrease in aggregate release flux and an increase in SR [Ca2+], 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 high-frequency/short-duration or low-frequency/long-duration Ca2+ sparks. The results are well correlated with recent experimental observations using pharmacological RyR inhibitors and clarify those aspects of the release-reuptake balance that are inherent to the coupling between local and global Ca2+ signals and those aspects that depend on molecular-level details. The model of Ca2+ sparks and homeostasis presented here can be a useful tool for understanding changes in cardiac Ca2+ release resulting from drugs, mutations, or acquired diseases.
- calcium signaling
- ryanodine receptor
- Markov chain
intracellular Ca2+ is a ubiquitous biological signal that serves diverse functions in many cell types. In individual cells, information can be conveyed by both “global,” or cell-wide, changes in [Ca2+] and by “local,” subcellular Ca2+ signals. Local signals are frequently caused by the release of Ca2+ 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,5-trisphosphate (IP3) receptors or ryanodine receptors (RyRs), and are observable experimentally as Ca2+ sparks (7) or puffs (51). When one or several of the channels in a release site are open, the [Ca2+] experienced by clustered channels is dramatically different from the [Ca2+] in the bulk cytosol.
In cardiac myocytes and many other cell types, local and global Ca2+ 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 [Ca2+] in the myoplasm or SR can influence the frequency, amplitude, and kinetics of local events. In ventricular myocytes, for instance, propagating Ca2+ waves emerge when spontaneous Ca2+ 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 [Ca2+]. An increase in RyR opening will cause a gradual decrease in SR [Ca2+] (46), whereas inhibition of RyR opening, over time, will lead to elevated SR [Ca2+] (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 Ca2+ 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 Ca2+ buffer protein that associates with and modulates RyRs (20, 33). Experiments in vitro have shown that CPVT-causing 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 Ca2+ responses can provide insights into disease pathophysiology and can potentially suggest novel therapies.
The difference in spatial scales between local and global Ca2+ signals, however, creates significant challenges for the development of mechanistic mathematical models. In particular, gating of RyRs depends on both myoplasmic and SR [Ca2+] (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 Ca2+ 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 Ca2+ sparks, but these have generally treated myoplasmic and bulk SR [Ca2+] as fixed boundary conditions. Conversely, modeling studies (4, 10, 40) focusing on cellular Ca2+ transients have usually used representations of SR Ca2+ release that do not account for the stochastic nature of the local events. Attempting to simulate Ca2+ 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 Ca2+ signals in permeabilized ventricular myocytes. The model accounts for the random generation and termination of spontaneous Ca2+ sparks, the resulting changes in myoplasmic and SR [Ca2+], 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 Ca2+ sparks followed by an increase in SR [Ca2+] and a partial recovery of spark frequency. Surprisingly, these authors found that prolonged exposure to tetracaine led to an increase in Ca2+ 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 [Ca2+]. More broadly, this model provides a powerful yet minimal framework for understanding how mutations, posttranslational modifications, or drugs can alter diastolic SR Ca2+ release in ventricular myocytes.
The minimal whole cell model of local and global Ca2+ responses developed here takes into account stochastic Ca2+-release site dynamics as well as the balance of release and reuptake fluxes leading to Ca2+ homeostasis in quiescent ventricular myocytes (Fig. 1). The model assumes that RyR Ca2+ 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 [Ca2+] (myoplasmic and SR), but these “domain” [Ca2+] are heterogenous throughout the cell, i.e., different release sites experience different domain [Ca2+]. Similar to previous work by Hinch and colleagues (16, 25, 26), we assume that when the number of open channels in a Ca2+-release site changes, the local [Ca2+] 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 Ca2+-release sites are coupled to the bulk myoplasmic and SR [Ca2+] in a manner that allows spontaneous Ca2+ sparks to change the balance of Ca2+ release, or “leak,” and reuptake by sarco(endo)plasmic reticulum Ca2+-ATPase (SERCA) pumps. The bulk myoplasmic and SR [Ca2+] determine the relationship between the number of open channels in a Ca2+-release site and the resulting domain [Ca2+], 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 Ca2+ handling is a novel aspect of our model formulation that has not been emphasized in previous work.
Ca2+-release site model.
In our model formulation, Ca2+-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) [Ca2+], denoted by cmyod, and closes with a Ca2+-independent transition rate: (1) where we assume cooperative Ca2+ binding. kco(cmyod)2 and koc are transition rates (in units of reciprocal time), and kco is an association rate constant (with units of concentration−2·time−1).
We assume that the local [Ca2+] experienced by the Ca2+-regulatory site of each channel depends only on the number of open channels at the Ca2+ release site [NO(t)] and the bulk Ca2+ concentrations (cmyo and csr), as described below. Because the channels are identical and indistinguishable, the state space for the N channel Ca2+-release site includes N + 1 states (0 ≤ NO ≤ N), as follows: (2) where the states (0, 1, … , N − 1, N) indicate the possible numbers of open channels and cmyod,n is the local myoplasmic domain [Ca2+] that applies when there are n open channels. The infinitesimal generator matrix, denoted by Q = (qij), that corresponds to Eq. 2 is tridiagonal: (3) where qij is the transition rate from state i to state j and the diamonds (◊) indicate a diagonal entry leading to a row sum of zero.
Once the local [Ca2+] (cmyod,n) that apply for any given number of open channels are specified, all of the statistical properties of the Ca2+-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: (4) where the row vector π = (π0, π1, … , πN), πi(t) is the probability of finding a release site in state i, and π(0) is the initial condition. The numerical solution of Eq. 4 can also be interpreted as providing the probability distribution for the state of a release site randomly sampled at time t from a large population that was initially prepared with a distribution of states given by π(0).
Myoplasmic and SR domain Ca2+.
As shown schematically in Fig. 1, the domain [Ca2+] for each release site (cmyod and csrd) are coupled to the bulk compartments via the myoplasmic and SR fluxes (Jmyo and Jsr) and coupled to one another through the release flux (Jrel). We assume these fluxes take the following form: (5) (6) (7) where cmyo and csr are the bulk myoplasmic and SR concentrations, γn = n/N indicates the fraction of channels at any given Ca2+-release site that are open, and vrel is the maximum release rate. Because the parameter vmyo is related to the exponential time constant for the decay of elevated myoplasmic domain [Ca2+] to the myoplasmic bulk [Ca2+] when an open release site closes, we will refer to vmyo as the rate of myoplasmic domain collapse (37). Similarly, the parameter vsr is the rate of luminal domain recovery that determines the exponential time constant for the relaxation of depleted junctional SR [Ca2+] as this compartment is refilled via Ca2+ translocation from the bulk SR (27).
Recall that Eqs. 2 and 3 require the specification of the local [Ca2+] that applies for any given number n of open channels (cmyod,n). Assuming the dynamics of domain Ca2+ are fast compared with the gating of RyRs, these domain [Ca2+] are found by balancing the fluxes Jreln = Jmyon and Jsrn = Jreln. Solving these equations simultaneously for cmyod,n and csrd,n yields the following: (8) (9) where (10)
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 (NO = n). Note that as the number of open channels increases, cmyod,n increases and csrd,n decreases, but always cmyod,n < csrd,n. The open circles in Fig. 2 show how an increase in the bulk SR [Ca2+] (csr) influences the relationship between the number of open channels in a Ca2+-release site and the resulting cytosolic (cmyod,n) and luminal (csrd,n) domain [Ca2+]. In particular, note that an increase in the bulk SR [Ca2+] (csr) leads to an increase in cytosolic domain [Ca2+] (cmyod,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 [Ca2+] are affected by the SERCA pump flux (Jpump), a passive leak from the SR to the myoplasm that is independent of release site activity (Jleak), as well as the fluxes from the myoplasmic domains to the bulk myoplasm (Jmyo) and from the bulk SR to the SR domains (Jsr). In a conventional Monte Carlo simulation involving a large number of Ca2+-release sites, concentration balance equations for the myoplasmic and SR [Ca2+] consistent with Fig. 1 are written as follows: (11) (12) where Jpm is plasma membrance flux, λsr = Vsr/Vmyo, Vmyo and Vsr are the effective myoplasmic and SR volumes (i.e., accounting for Ca2+-buffering capacity), and JmyoT and JsrT are total fluxes obtained summing over all release sites. Under the assumption of fast domain Ca2+, these fluxes can be expressed as a sum over all possible release site states as follows: (13) (14) where fn(t) are random variables denoting the fraction of release sites with n open channels, 0 ≤ fn ≤ 1, Σn = 0N fn = 1, cmyod,n and csrd,n are given by Eqs. 8–10, and the rate constants vmyoT and vsrT are proportional to vmyo and vsr but scaled by the total number of release sites. In a conventional Monte Carlo approach to simulation of a whole cell model of local and global Ca2+ dynamics, the ODEs for bulk myoplasmic and SR Ca2+ are integrated while bidirectionally coupled (through fn in Eqs. 11–14) to a stochastic simulation of a large finite number of Markov chains (Eq. 2), each one of which represents a Ca2+-release site.
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 fn in Eqs. 13 and 14. Unfortunately, this approach is invalid because the assumption of rapid Ca2+ 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 Ca2+ (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) [Ca2+]: (15) (16) where the total release flux JrelT is given by the following: (17) where γn = n/N, cmyod,n and csrd,n are given by Eqs. 8–10, πn is the probability that a randomly sampled release site has n open channels, and, as mentioned above, π = (π0, π1, … , πN) is found by integrating Eq. 4. The total myoplasmic (ĉmyo) and SR (ĉsr) [Ca2+] that solve Eqs. 15 and 16 are the following sums of the bulk and domain concentrations weighted by effective volume ratios: (18) (19) In these definitions, c̄myod and c̄srd are the average myoplasmic and SR domain [Ca2+] that would be obtained upon randomly sampling release sites from within the cell, that is: (20) (21) The effective volume ratios that appear in Eqs. 18 and 19 are given by Λmyod = Vmyod,T/Vmyo, Λsrd = Vsrd,T/Vmyo, and λsr = Vsr/Vmyo, where Vmyo and Vsr are the effective myoplasmic and SR volumes and Vmyod,T and Vsrd,T are the effective volumes of the aggregated myoplasmic and SR domains, respectively.
Because the experimental observations that are of primary relevance to this report involve permeabilized cells (52), the plasma membrane flux (Jpm) was chosen to be as follows: (22) where kpm was chosen to be large enough to “clamp” the bulk myoplasmic [Ca2+] (cmyo) at the level of the extracellular bath (cext = 0.1 μM in the standard parameter set). Even so, the total myoplasmic [Ca2+] (ĉmyo, Eq. 18) that solves Eq. 15 is not fixed, because this concentration includes Ca2+ in the myoplasmic domains. The fluxes between the bulk myoplasm and bulk SR that occur in Eqs. 15 and 16 include Ca2+ reuptake by SERCA pumps (Jpump): (23) and a passive leakage flux (Jleak): (24) This minimal model of the relationship between Ca2+ sparks and Ca2+ homeostasis was implemented in Matlab (The Mathworks) running on a 1.67-GHz Power PC with 1-GB memory. The model ODEs are stiff and were integrated using Matlab's built-in function ode15s using an adaptive time step and relative and absolute tolerances of 10−3 and 10−6.
Summary and significance of the model.
Although minimal in nature, the whole cell model of local and global Ca2+ signaling that is the focus of this report accounts for the changes in myoplasmic and SR [Ca2+] mediated by the balance of stochastic release and reuptake by the SR and the feedback of myoplasmic and SR [Ca2+] on spark frequency. As discussed in the Introduction, previously published models of Ca2+ signaling in cardiac myocytes that include the stochastic release of SR Ca2+ either have not included bidirectional coupling between local Ca2+ release and global Ca2+ 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 Ca2+, and the resulting balance of bulk [Ca2+] in permeabilized ventricular myocytes. It is also the first model of Ca2+ sparks and homeostasis that bypasses Monte Carlo simulation by assuming both a large number of Ca2+-release sites and rapid Ca2+ domain dynamics, resulting in a minimal formulation that facilitates parameter studies.
The minimal whole cell model of local and global Ca2+ 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) [Ca2+] (Eqs. 15 and 16). The additional N + 1 ODEs (Eq. 4) account for the dynamics of a large number of Ca2+-release sites, each composed of N two-state 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 (cmyod,n) and SR (csrd,n) domain [Ca2+] and the number of open channels (Eqs. 8–10). Note that the fluxes JrelT, Jpump, and Jleak are functions of cmyo and csr, 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 Ca2+-release sites, we do not have to specify a precise number. To see this, note that the domain concentrations cmyod,n and csrd,n do not depend on the number of release sites in the cell (M) when the rate constants are defined by vrel = vrelT/M, vmyo = vmyoT/M, and vsr = vsrT/M (Eqs. 8 and 9). Because the rate vrelT that appears in Eq. 17 does not correspond to release through one release site but rather the entire population, it is convenient to specify cmyod,n and csrd,n using Eqs. 8–10 with the replacement of vrelT, vmyoT, and vsrT for vrel, vmyo, and vsr.
The minimal model of local and global Ca2+ signaling includes 14 parameters, far fewer than most mathematical models of Ca2+ handling in cardiac myocytes (see the Supplemental Table in the Supplemental Material). Some parameters [such as the effective volume ratios λsr, Λsrd, and Λmyod and the SERCA pump maximum rate (vpump) and dissociation constant (kpump)] 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 kpm is unimportant so long as there is rapid equilibration of bulk myoplasmic Ca2+ with the extracellular [Ca2+] (cext). 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 Ca2+ 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 three-dimensional electron tomography (2, 23). The most important of the model parameters [the kinetic parameters for the stochastic gating of the two-state RyR (kco and koc) and the rate constants for Ca2+ release vrelT, myoplasmic domain collapse vmyoT, and luminal domain recovery vsrT] 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 [Ca2+] (cmyo), the average duration of a spontaneous Ca2+-release event is on the order of 20 ms, similar to the observed rise time of Ca2+ sparks (6, 7). The Ca2+ spark rate with cmyo = 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 two-dimensional 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 cross-sectional area of 100 × 30 μm. Consistent with experiment, an increase in myoplasmic [Ca2+] in the model leads to an increase in the spontaneous Ca2+ spark rate.
RyR open probability and spontaneous Ca2+ sparks.
The minimal model of local and global Ca2+ signaling that is the focus of this report simulates stochastic Ca2+ release by clusters of RyRs and the resulting whole cell Ca2+ 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 Ca2+ 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 Ca2+-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 two-state RyR model is given by τC = 1/kco(cmyod)2 (Eq. 1), we simulated the application of tetracaine to permeabilized ventricular myocytes by decreasing the rate constant kco, which influences the stochastic dynamics of the Ca2+-release sites (Eq. 2) in the minimal whole cell model. We wanted to understand how the simulated application of tetracaine influences the dynamics of Ca2+ 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 Ca2+ activation rate constant kco. The circles show the result of two particular simulations: one corresponding to the standard parameter values (kco = 4.5 μM−2·s−1; solid circles) and the other corresponding to the simulated application of tetracaine (kco = 0.5 μM−2·s−1; open circles). In the latter case, the value of kco was chosen so that the single channel probablity (Popen), given by the following: (25) decreased by 88% upon the action of tetracaine, consistent with experiments in which 0.7 mM tetracaine was applied (52).
Figure 3A shows that the simulated application of tetracaine led to increased bulk SR [Ca2+] in the whole cell model (csr = 342 to 1,112 μM; compare solid and open circles). As expected, steady-state bulk SR [Ca2+] 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 [Ca2+], maximum bulk SR [Ca2+] (csr) asymptotically approaches 2.5 mM when association rate constant kco is very small. When a nonspecific passive leak was not included in the model, csr 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 (fO) 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 [Ca2+] and the interaction between RyRs combine to attenuate the decrease in channel activity occurring during the simulated application of tetracaine. That is, increased SR [Ca2+] increases the driving force during stochastic Ca2+-release events and elevates myoplasmic domain [Ca2+] (compare solid and open circles in Fig. 2B). The fraction of open channels can be calculated as fO = E[NO]/N, where E[NO] is the expected value of the number of open channels in a randomly sampled release site: (26) Figure 3, C and D, shows the frequency and duration of spontaneous Ca2+ sparks occurring in the whole cell model as a function of the parameter kco. The presence or absence of Ca2+ sparks is assessed by calculating the following Ca2+ spark score: (27) where (28) The Ca2+ spark score takes values from 0 to 1, and a score of 0.2 or higher indicates robust sparks (38). The solid lines in Fig. 3, C and D, show that Ca2+ sparks were observed when the RyR Ca2+ activation rate constant (kco) was between 0.04 and 322 μM−2·s−1, a range spanning four orders of magnitude.
With spark initiation defined as a release site reaching a threshold number of open channels (NO = 4 to 5 transition) and spark termination defined as all channels closing (NO = 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 Ca2+ spark frequency but increased mean Ca2+ spark duration, consistent with experimental observations (52). While bulk SR [Ca2+] (csr) 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 Ca2+-release events exhibited by a Ca2+-release site in the standard whole cell simulation. Spontaneous Ca2+ sparks were simulated using Gillespie's method (14). The parameters used correspond to the solid circles shown in Fig. 3 (kco = 4.5 μM−2·s−1) and resulted in robust Ca2+ sparks (score = 0.51) for the bulk concentrations (cmyo and csr) of the equilibrated whole cell model. Note the high frequency of spontaneous release events, including five sparks (NO ≥ 5) and a large number of smaller release events. These small release events, termed “Ca2+ quarks” (32), would not be detectable with standard confocal microscopy and would therefore contribute to “invisible” SR Ca2+ 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 Ca2+-release events during the simulated addition of tetracaine. The parameters used correspond to the open circles shown in Fig. 3 (kco = 0.5 μM−2·s−1). While robust Ca2+ 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 long-duration 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 [Ca2+], Fig. 4C shows a control simulation using the single channel RyR parameters shown in Fig. 4B (kco = 0.5 μM−2·s−1) with bulk SR [Ca2+] “clamped” at the value shown in Fig. 4A (csr = 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 [Ca2+] 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 Ca2+ 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 Ca2+ release due to spontaneous sparks.
As discussed above, the SR [Ca2+] overload induced by tetracaine led to higher myoplasmic domain [Ca2+] (compare open and solid circles in Fig. 2A), higher SR domain [Ca2+] (Fig. 2B), and higher release flux (JrelT, Eq. 17) for any given number of open channels. The resulting changes in the dynamics of Ca2+-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 [Ca2+] is monotone decreasing (Fig. 3A).
During experimental observations of spontaneous Ca2+ release, small-amplitude 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 Ca2+ activation rate constant kco. Figure 5D shows aggregate flux JrelT 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 NO = 1 as well as NO = 6 or 7. Tetracaine suppressed the proportion of Ca2+ released through sites with seven or fewer open channels (NO ≤ 7), whereas release mediated by sites with eight or more open channels (NO ≥ 8) increased slightly in the presence of tetracaine. These observations are consistent with the increased probability of long-duration Ca2+ 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 Ca2+ 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 NO = 0 → 1 transition, whereas spark amplitude was the integrated stochastic release flux (Eq. 5) before spark termination via a NO = 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” Ca2+-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 Ca2+ 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 Ca2+ 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 steady-state dynamics of the whole cell model, Fig. 7 shows transient effects upon bulk SR [Ca2+] (csr), 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 Ca2+ release caused a slow increase in SR load that ultimately increased the mean spark duration to 26.8 ms, consistent with the steady-state 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 [Ca2+] 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 Ca2+ homeostasis.
Figure 8 shows a summary of 2,500 calculations of the stationary dynamics of the whole cell model as a function of vsrT and vmyoT. These parameters control the rate of “diffusion” or translocation of Ca2+ 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 Ca2+-mediated coupling between RyRs and the extent of SR domain depletion during sparks (recall Eqs. 8–10). Because vsrT and vmyoT 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 steady-state bulk SR [Ca2+] and mean spark frequency and duration.
Figure 8A shows that bulk SR [Ca2+] is an increasing function of rate of myoplasmic domain collapse (vmyoT) and a decreasing function of the rate of SR domain recovery (vsrT). When vmyoT is very large or vsrT 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 [Ca2+]. Figure 8, B and D, shows that the fraction of open channels and mean spark duration were much more sensitive to vmyoT than to vsrT. When vmyoT is small, sparks are extremely long because subspace [Ca2+] remains elevated after RyRs close (cf. Refs. 27 and 37). Fast SR refilling alone is not sufficient to induce long-duration Ca2+ sparks of the sort observed upon the simulated application of tetracaine (cf. Fig. 4B).
Robust Ca2+ 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 Ca2+ spark behavior is a consequence of the homeostatic changes in bulk SR [Ca2+] accounted for in our model formulation; when these simulations were repeated with bulk myoplasmic and SR Ca2+ concentrations fixed at baseline values (cmyo = 0.1 μM and csr = 342 μM), the range of domain rate constants leading to Ca2+ sparks was considerably smaller (Fig. 9B).
Figure 10 shows stationary dynamics of the whole cell model as a function of bulk myoplasmic [Ca2+] (cmyo) and maximum rate of the SERCA pump (vpump), two parameters that can be easily manipulated in experiments. Bulk SR [Ca2+] was (as expected) an increasing function of vpump; however, the SR Ca2+ load was a biphasic function of cmyo (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 [Ca2+] (Fig. 10C). Spark duration was a biphasic function of cmyo (Fig. 10D), consistent with the biphasic effect of this parameter on SR load (Fig. 10A). Bulk myoplasmic [Ca2+] (cmyo) significantly changed the functional dependence of the stationary dynamics of the whole cell model on the RyR Ca2+ activation rate constant (kco). On the other hand, when the application of tetracaine was simulated by comparing kco values that correspond to a decrease in RyR activity from fO = 9.6 × 10−4 to 2.0 × 10−4 (Table 1), we observed increased SR load, decreased spark frequency, and increased spark duration for cmyo 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 Ca2+ activation rate constant kco, which reduces the open probability (Eq. 25) of the RyR model (Eq. 1) by increasing the mean closed dwell time, τC = 1/kco(cmyod)2. However, the open probability of the RyR can also be reduced by increasing the rate constant koc, thereby decreasing the mean open dwell time (τO = 1/koc). 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 Ca2+ spark frequency and mean spark duration as a function of the RyR kinetic constants kco and koc 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 = koc/kco, 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 [Ca2+] also depend only on RyR open probability (Fig. 11A). SR [Ca2+] 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 kco, koc, or both (diamonds). A reduction in kco (increase in τC), analogous to effect of tetracaine, led to fewer sparks but an increase in spark durations. Conversely, an increase in koc (decrease in τO), analogous to the effect of flecainide, caused more sparks with decreased mean duration. Examples of Monte-Carlo simulations under these different conditions are shown in Fig. 12. Similarly, we can assume that low-dose caffeine leads to an increase kco accompanied by a smaller decrease in koc (+ symbol) (29). The model predicted that this would cause a decrease in SR [Ca2+], 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 Ca2+ 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 steady-state bulk SR [Ca2+], are less dependent on the kinetic details of the perturbation but remain determined by equilibrium quantities such as the dissociation constant for Ca2+ binding to the RyR.
This report presents a minimal whole cell model of local and global Ca2+ signals in quiescent ventricular myocytes. The modeling formalism accounts for the effect of random spontaneous Ca2+ sparks, changes in bulk myoplasmic and SR [Ca2+] mediated by the balance of stochastic release and reuptake by the SR, and feedback of myoplasmic and SR [Ca2+] 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 Ca2+ are similar to assumptions made in previously published local control models that represent the stochastic dynamics of a large number of Ca2+-release sites (16, 25), but these previous studies have not made a distinction between domain and bulk SR Ca2+ 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 Ca2+ release, and the resulting balance of bulk myoplasmic and SR [Ca2+] in quiescent ventricular myocytes. Because of the computational challenge of large-scale 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 Ca2+-release sites and rapid Ca2+ domain dynamics, resulting in a minimal formulation that facilitates parameter studies.
The minimal model of local and global Ca2+ 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 Ca2+ sparks that is followed by an increase in bulk SR [Ca2+], partial recovery of spark frequency, and an increase in Ca2+ 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 steady-state SR [Ca2+], 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, Ca2+ spark characteristics, and balance of bulk myoplasmic and SR [Ca2+] (Figs. 3⇑–5). For example, simulations showing that mean spark duration is a biphasic function of the RyR Ca2+ activation rate constant kco 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 Ca2+-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 Ca2+ balance is largely determined by RyR open probability (a function of K = koc/kco). Conversely, the frequency and duration of Ca2+ sparks are sensitive to RyR open and closed dwell times, which are determined by kco and koc independently (see Fig. 11).
The Ca2+-release site model used here is quite simple and assumes “instantaneous mean field coupling” (38) of N = 10 two-state 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 Ca2+-release site model is not so large that integrating Eq. 4 is impractical. Our use of a two-state RyR in this first study of the relationship between Ca2+ 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 Ca2+ 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 Ca2+-dependent inactivation and/or sensitization by SR [Ca2+]. The effect of these well-established aspects of RyR Ca2+ regulation is beyond the scope of this work.
The Ca2+ activation process in the RyR model is mediated by myoplasmic domain Ca2+ (cmyod, Eq. 9), and it would be possible to augment the model to include the sensitization of RyRs by SR domain Ca2+ (csrd, Eq. 9). However, the assumption of instantaneous coupling of RyRs (i.e., rapid equilibration of myoplasmic and SR domain Ca2+) may not work well when luminal Ca2+ plays a role in spark termination (27, 42). A computational study of the contribution of luminal RyR regulation to the bidirectional coupling of spontaneous Ca2+ release and Ca2+ homeostasis will likely require more subtle mathematical formulations, similar to the probability density and moment closure techniques that accelerate “local control” simulations of high-gain graded Ca2+ release in voltage-clamped cardiac myocytes (49, 50). These representations of heterogenous domain Ca2+ concentrations associated with a large number of Ca2+-release sites remain valid even when the dynamics of SR domain Ca2+ 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 Ca2+ domains rapidly equilibrate with myoplasmic domain and bulk SR Ca2+, in which case there is a negligible variance of SR domain [Ca2+] for any given Ca2+-release site state.
In preliminary work, we studied whole cell responses using RyR models that included both Ca2+ activation and inactivation by myoplasmic domain Ca2+, for example, the following three-state model: (29) which includes a long-lived closed state (R). As expected, release sites composed of such channels exhibit Ca2+ sparks and homeostasis similar to that observed for the two-state RyR when the dissociation constant Kinact = kro/kor is large (see supplemental Fig. S3). While fast Ca2+ inactivation (large kro and kor with fixed Kinact) is associated with decreased spark duration and increased spark frequency, extremely fast Ca2+ inactivation can preclude sparks (supplemental Fig. S3, C and D). Provided Ca2+ inactivation is sufficiently fast, decreasing Kinact leads to increased SR Ca2+ load, decreased RyR activity, increased spark frequency, and decreased spark duration (supplemental Fig. S3, A–D). Interestingly, decreasing the Ca2+ activation rate constant (kco) in the three-state model with Ca2+ inactivation (Eq. 29) to simulate the application of tetracaine may lead to longer or shorter duration sparks depending on whether the rate of Ca2+ inactivation is slow or fast, respectively, despite the fact that neither RyR activity nor SR load are strongly affected by the rate of Ca2+ inactivation in either condition (see supplemental Fig. S4). This sensitivity of the stochastic dynamics of Ca2+ release to the rate of Ca2+ inactivation is consistent with results obtained using Ca2+-release site models that do not account for Ca2+ homeostasis (18). However, the significance of these observations is unclear given recent experimental results showing that, at even 50 μM myoplasmic [Ca2+], inactivation is unable to suppress SR Ca2+ release in permeabilized myocytes (45).
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 (vsrT). Consistent with the hypothesis of Zima et al. (52), an increase in vsrT from 10 to 50 s−1 caused a slight increase in baseline Ca2+ 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 Ca2+ 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 Ca2+ sparks is considerably smaller when bulk myoplasmic and SR [Ca2+] are fixed at baseline values (compare Fig. 9, A and B), we conclude that this robust Ca2+ spark behavior is a consequence of homeostatic changes in bulk SR [Ca2+] 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 Ca2+ sparks and discourage the alternatives: SR leak mediated by quarks or tonically active release sites. These observations underscore the importance of accounting for global Ca2+ balance in models of localized Ca2+-release events.
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.
No conflicts of interest, financial or otherwise, are declared by the author(s).
Some of the results of this study have previously appeared in abstract form (22).
↵1 Supplemental Material for this article is available online at the American Journal of Physiology-Heart and Circulatory Physiology website.
- Copyright © 2010 the American Physiological Society