Calcium sparks in cardiac myocytes are brief, localized calcium releases from the sarcoplasmic reticulum (SR) believed to be caused by locally regenerative calcium-induced calcium release (CICR) via couplons, clusters of ryanodine receptors (RyRs). How such regeneration is terminated is uncertain. We performed numerical simulations of an idealized stochastic model of spark production, assuming a RyR gating scheme with only two states (open and closed). Local depletion of calcium in the SR was inevitable during a spark, and this could terminate sparks by interrupting CICR, with or without assumed modulation of RyR gating by SR lumenal calcium. Spark termination by local SR depletion was not robust: under some conditions, sparks could be greatly and variably prolonged, terminating by stochastic attrition–a phenomenon we dub “spark metastability.” Spark fluorescence rise time was not a good surrogate for the duration of calcium release. Using a highly simplified, deterministic model of the dynamics of a couplon, we show that spark metastability depends on the kinetic relationship of RyR gating and junctional SR refilling rates. The conditions for spark metastability resemble those produced by known mutations of RyR2 and CASQ2 that cause life-threatening triggered arrhythmias, and spark metastability may be mitigated by altering the kinetics of the RyR in a manner similar to the effects of drugs known to prevent those arrhythmias. The model was unable to explain the distributions of spark amplitudes and rise times seen in chemically skinned cat atrial myocytes, suggesting that such sparks may be more complex events involving heterogeneity of couplons or local propagation among sub-clusters of RyRs.
Resting, mammalian cardiac muscle produces local calcium release events known as calcium sparks (Cheng et al., 1993). It is widely accepted that these events are caused by coordinated CICR (Fabiato, 1983) among local clusters (“couplons”; Stern et al., 1997, 1999) of calcium release channels (RyRs), usually located in a dyadic cleft between the SR and the sarcolemma. This creates a paradox. If spark release is caused by local, regenerative positive feedback of CICR, what causes or allows it to terminate? Early on, it was believed that a time- and Ca2+-dependent inactivation of RyRs was involved (Stern and Cheng, 2004), but studies of isolated RyRs of the cardiac isoform in lipid bilayers have not demonstrated a sufficiently powerful inactivation process (although it remains possible that such inactivation exists in vivo or is mediated by allosteric interactions, as suggested by Stern et al., 1999). It is now generally accepted that local depletion of calcium in the terminal cisternae of the junctional SR (JSR) is involved in spark termination. A common hypothesis (Lukyanenko et al., 1998; Terentyev et al., 2002) is that when [Ca2+]JSR falls below a critical threshold (often suggested to be ∼50% of the resting level), sparks terminate because of deactivation of RyRs, either by direct action of a lumenal calcium–sensing site on the RyR or mediated by the calcium-buffering protein calsequestrin (Casq) located in the JSR.
In this paper, we examine this hypothesis semi-quantitatively by means of numerical modeling. Specifically, we generate statistical distributions of spark properties from a minimal model of stochastic CICR in spatially extended RyR clusters, and try to determine whether, and under what conditions, observed spark statistics could be explained by such a model. We find that spark termination can be produced by local JSR depletion, with or without the hypothesized regulation of RyR gating by lumenal calcium. However, the robustness of spark termination depends on kinetic relationships among rates of RyR gating and of SR calcium recycling. In other words, the ability of depletion to terminate sparks is not absolute but depends on the kinetics of RyR gating and JSR refilling. An important consequence of this result is prediction of a novel regimen of sparks, dubbed “metastable sparks,” which persist for a prolonged period before terminating by stochastic attrition. Anything that makes the JSR refill more quickly or causes a greater sensitivity of the RyR to calcium, or a slower RyR closure rate (longer open dwell time), will favor the metastable spark regimen.
MATERIALS AND METHODS
Calcium sparks were simulated by an extension of a stochastic algorithm described previously (Stern et al., 1999). In brief, the algorithm generated a statistical ensemble of 10,000 realizations of a model couplon consisting in a square array of RyRs located on a 30-nm lattice in a dyadic cleft. Sparks were initiated by randomly opening one RyR in the couplon. Dyadic cleft [Ca2+] was calculated by a reaction–diffusion equation discretized on a two-dimensional 10-nm grid as described previously (Stern et al., 1999). In addition, the differential system now included a JSR compartment, with instantaneous buffering by Casq, and, optionally, a chain of free SR (FSR) compartments without buffering. The JSR compartment (Fig. 1) was coupled to an infinite pool of FSR calcium through either a single-diffusion resistance (Fig. 1 A) or a chain of 100 compartments simulating the local FSR of the sarcomere (Fig. 1 B). The joint gating of the RyRs in a couplon was simulated by a Monte Carlo algorithm that produces an exact realization of the variable-rate, continuous-time Markov process driven by the evolving cleft calcium distribution (Stern et al., 1997).
For this study, we confined ourselves to a minimal gating scheme in which a single RyR has only two states, closed and open, with no inactivation. The transition rates were given as a function of local cleft and lumenal calcium concentrations by:(1.1)(1.2)where hry is the cooperativity of RyR opening by calcium, ko is the opening rate constant of the RyR (dimension mM−hry ms−1), and the parameter sloc0 determines the relative change in opening rate as lumenal calcium is decreased from resting. The usual value of sloc0, 0.1 mM, implies an 11-fold decrease of opening rate as SR calcium falls from 1 mM to zero. For the simulations presented here, hry, the cooperativity of RyR activation by cytosolic calcium, was taken to be 2, in keeping with many studies of isolated RyRs in lipid bilayer (see Qin et al., 2009, for a graphical summary of many results). We did not include allosteric coupling (“coupled gating”; Marx et al., 2001) among adjacent RyRs, as demonstrations of this phenomenon in RyR2 have been few and controversial.
To compare simulations with experimental spark fluorescence data, the total RyR flux from each couplon was taken as source for a spherically symmetrical diffusion and buffering computation, and convolved with an empirical confocal point–spread function, as described previously (Santiago et al., 2010). Histograms of the distribution of observable in-focus spark properties (e.g., rise time, amplitude) were constructed from the ensemble of detected sparks out of 10,000 individual spark simulations generated in parallel. Sparks were considered to be “detected” if they satisfied the criteria that peak F/F0 > 1.2, with a peak rate of rise >0.2 ms−1.
To examine further the mechanisms of spark initiation and termination, we computed the behavior of a simple, deterministic two-compartment continuum model of a couplon (“toy model”). Further details of the model structures, parameters, and implementation are given in the Appendix.
Online supplemental material
The online supplemental material includes a description of structure and implementation of the stochastic model, RyR gating estimation of Dup, and deterministic simulations of sparks/blinks. Fig. S1 shows simulated central spark fluorescence using a variety of different endogenous buffer parameters taken from published papers, assuming the same, decaying release flux through a fixed number of open RyRs. The calculated spark amplitude varies widely depending on the assumptions. Fig. S2 shows JSR calcium computed from the full, deterministic simulation using a constant square-wave release flux, compared with a single-compartment model using the value of Dup estimated from steady-state diffusion as described in the supplemental text. Fig. S3 shows that peak fluorescence is reached early, regardless of the length of RyR opening. Even for large rates of JSR refilling, depletion causes fluorescence to decline. Table S1 lists typical model parameters. Videos 1 and 2 show the RyR gating dynamics of a normal spark and a metastable spark, respectively. The online supplemental material is available at http://www.jgp.org/cgi/content/full/jgp.201311034/DC1.
How big is a spark?
To determine whether it is feasible to compare simulated spark amplitudes with experimental ones, we performed a series of deterministic spark fluorescence computations, each using the same calcium source, but assuming four different systems of cytosolic buffers, each of which had been used in a published spark modeling or analysis paper (Smith et al., 1998; Soeller and Cannell, 2002; Santiago et al., 2010; Shkryl et al., 2012). As shown in Fig. S1, the peak amplitude of the simulated sparks ranged over an order of magnitude depending on which buffering system was assumed. We conclude that no absolute comparison can be made between the amplitudes of simulated and observed sparks unless the parameters of cytosolic diffusion and buffering have been measured in the same preparation from which the sparks were recorded. Henceforth, we concentrate on statistical properties of spark populations that are independent of absolute amplitude.
Estimation of JSR refilling
Modeling sparks requires estimation of the rate at which calcium from the FSR network refills the JSR release terminal. The detailed architecture and topology of the cardiac FSR are not well known. The simplest model assumes that the JSR refills from an infinite FSR pool through a single effective diffusion resistance. This resistance–scaled to JSR volume–may be characterized by the parameter Dup (dimension ms−1), which is the rate constant of JSR refilling that would pertain in the absence of Casq. Estimation of Dup from first principles turns out to be problematic.
A simple model is to treat the FSR as a fine-grained random network occupying a volume fraction fvfsr (∼3.5%; Shannon et al., 2004) of the cell. By a theorem of stereology, this is also the fraction of any surface crossed by the FSR lumen. As a result, the diffusion equation for FSR calcium is formally the same as that in the cytosol, except that boundary fluxes are reduced by the factor fvfsr. If we assume spherical geometry, with the JSR occupying a sphere of radius rjsr and volume Vjsr, we can estimate Dup as the steady-state diffusional conductance (flux over concentration difference) between rjsr and infinity, divided by the JSR volume. Expressing rjsr in terms of the fractional volume fvjsr of JSR and the length and radius of the sarcomere gives (see supplemental text):(1.3)where Dsr is the calcium diffusion coefficient in the SR lumen, estimated to be 0.06 µm2/ms (Picht et al., 2011).
If we assume that the JSR occupies 0.35% of the volume of a cylindrical sarcomere of length 2 µm and diameter 0.8 µm, this gives Dup = 0.70765 ms−1. Although this calculation was based on a steady-state diffusion gradient, Fig. S2 shows that a single-compartment JSR simulation using this value, including 30-mM Casq-binding sites, reproduces well the JSR depletion and refilling of a full, deterministic spark/blink simulation. The problem is that, in both simulations, the refilling phase (∼10-ms time constant in the presence of Casq) is still much too rapid when compared with experimentally observed blinks, which have time constants of ∼200 ms in rabbit (Picht et al., 2011), or experimental spark restitution, which recovers approximately as a 90-ms exponential in rat (Sobie et al., 2005; Ramay et al., 2011). Alternatively, we may estimate Dup from these observed time constants, by assuming that the observed time constant τ is the exponential time constant of a one-compartment JSR with Casq, linearized about the resting state (see supplemental text):(1.4)This gives values for Dup in the range of 0.05–0.15 ms−1, considerably less than that estimated above from the random network model.
A more detailed simulation of FSR diffusion has been given by Picht et al. (2011), in which the FSR network was modeled as linear segments connecting JSR terminals in a three-dimensional grid. We ran their model with an assumed 20-ms square-wave release current. Although the resulting JSR calcium (not depicted) had a nonexponential time course, the time to e-fold recovery was ∼30 ms in the presence of Casq, corresponding to Dup ≈ 0.3 ms−1. Picht et al. (2011) hypothesized that the much slower observed blink recovery was caused by continuing RyR release on the tail of a spark. As we show below, this interpretation is problematic.
In summary, there is a discrepancy between the observed rates of blink refilling and spark restitution versus simple models of intra-lumenal diffusion in the SR. Presumably, the discrepancy will be accounted for either by a more accurate measurement of the intra-SR calcium diffusion coefficient or by further details of the nano-anatomy of the SR, for example, a constriction or bottleneck in the connection of the JSR to the FSR (Brochet et al., 2005). In the meantime, we will assume that Dup lies somewhere in the range of 0.05–1 and consider the consequences of these possible values for spark termination.
Lumenal depletion is unavoidable
Most experimental studies have estimated that only ∼50% of JSR calcium is released during a spark. In the Appendix we calculate how much lumenal [Ca2+] would be reduced, in steady state, by an indefinitely prolonged opening of nopen RyRs. Fig. 2 shows that the steady-state JSR depletion produced by even a single, open RyR would be >50% even at the highest Dup, corresponding to a refilling time of ∼10 ms, much faster than observed. Obviously, the large number of RyRs that participate in a spark would be expected to produce much more severe depletion, especially with a more realistic refilling rate. Therefore, major JSR depletion is inevitable unless a spark were terminated very early by a very powerful inactivation process. This suggests that the degree of local JSR depletion (“blink”) during a spark may have been underestimated by confocal microscopy.
The meaning of “spark duration”
It is often stated that the rise time of spark fluorescence at the center of the spark is an approximation of the duration of calcium release during the spark. This is not true for the rapidly declining release flux from a depleting JSR terminal. As shown in the supplemental text, rapid depletion of the JSR causes the amplitude and timing of the fluorescence peak to be controlled by the amount of calcium in the JSR, independent of the duration of the release. From a mechanistic point of view, it would be better to define spark duration as the time elapsed between the opening of the first RyR that triggers the spark and the closing of the last RyR causing extinction of CICR. Fig. 3 A shows the time courses of the number of open RyRs; JSR free calcium and dye fluorescence for a representative, stochastic, and simulated spark; as well as histograms of fluorescence rise time and spark extinction time from “detected” sparks out of 10,000 stochastic sparks attempted. There is clearly no similarity between these two measures of spark duration, and, as shown by the scatter plot (Fig. 3 B), there is very little relationship between them. Unfortunately, spark extinction time is not reliably observable. Once JSR depletion has occurred, the tail of calcium release flux is likely to be lost in the noise inherent in inverting the spark fluorescence signal.
Implicit in the histogram of Fig. 3 is that sparks do self-extinguish in a reasonable time. In other words, under the assumed conditions, JSR depletion provides a sufficient spark termination mechanism without requiring any RyR inactivation process. As shown in Fig. 4, this termination does not require the dependence of RyR gating on lumenal [Ca2+]. Spark extinction time was only modestly increased when lumenal activation was turned off. This shows that the interruption of positive feedback of CICR among RyRs within a couplon–by reduction in RyR flux as the JSR becomes locally depleted–is a sufficient mechanism to terminate sparks (as also shown by two models recently published: Gillespie and Fill, 2013, and Laver et al., 2013).
The robustness of spark termination, however, is strongly affected by the rate of refilling of the JSR from the neighboring FSR. The sparks in Figs. 3 and 4 were generated under conditions when Dup was set to 0.15 ms−1, as inferred from empirical restitution and blink relaxation rates by Eq. 1.4. If Dup is increased to 0.7 ms−1–the value predicted for a random network model of the SR–then a substantial fraction of sparks fails to terminate early, and continues to “smolder” for a prolonged period, forming a straggling tail on the histogram of extinction times (Fig. 5). We refer to this predicted behavior as “spark metastability.”
An obvious (but insufficient) explanation for the failure of these sparks to terminate would be that continuous, rapid refilling of the JSR prevents [Ca2+]JSR from being reduced to the “threshold” for spark termination. As Fig. 2 shows, even a modest number of open RyRs can severely deplete JSR calcium, even at the highest refilling rate; as shown below, this occurs rapidly on a time-scale short compared with spark duration. Another explanation is needed.
A toy model of the couplon
To examine the spark extinction mechanism more closely, we constructed a deterministic “toy” model of a single couplon. We treat the couplon as a “bag” of RyRs sharing a single, common cleft compartment whose calcium is in instantaneous steady state. The number of open RyRs at a given instant is treated as a continuous variable. The flux through an open RyR is proportional to the difference between [Ca2+]JSR and [Ca2+]cleft. [Ca2+]cleft is determined by the balance between aggregate RyR flux into the cleft and diffusion out of the edges of the cleft: we estimate [Ca2+]cleft as the average value of [Ca2+] in a circular, thin (two-dimensional) cleft if the total RyR flux were deposited at the center. [Ca2+]cleft then depends only on the flux, the calcium diffusion coefficient, and the height (but not the diameter) of the cleft (see Appendix). The instantaneous “number” of open RyRs varies at a rate equal to the difference between the opening rate and closing rate as determined by the kinetics of the RyR gating scheme, Eqs. 1.1 and 1.2.
With these simple assumptions, the couplon becomes a deterministic dynamical system with two state variables, nopen and casr, governed by a pair of differential equations that appear rather daunting once all the simple pieces above have been assembled (see Appendix):(1.5)(1.6)
We consider first Eq. 1.6 for nopen, with casr held constant. Taking the RyR calcium cooperativity hry to be 2, this equation has three fixed points (i.e., equilibrium values of nopen at which its rate of change vanishes). One stable fixed point describes a resting condition with nopen ≈ 0, a tiny value because of openings stimulated by resting background calcium. Another stable fixed point describes the condition where nopen ≈ nry and most RyRs are open, kept open by the large cleft [Ca2+] caused by their currents. Between these two lies an unstable fixed point, describing a threshold (at a given casr) for regeneration of CICR. As shown in Fig. 6 A, if nopen is started below the threshold (green curve), it will fall to the closed fixed point, whereas if it is above (red), it will increase (initially exponentially) to converge to the fully open state. The value of the threshold depends inversely on casr, as shown in Fig. 6 B.
We can now describe the sequence of events when Eqs. 1.5 and 1.6 are coupled. Initially, the JSR is fully loaded, so the threshold for CICR is very low (Fig. 6 B), reflecting the fact that RyRs are highly sensitive and RyR current is high, favoring regeneration of CICR. If nopen is raised above this threshold (shown in Fig. 6 A as a “bifurcation point” at time 0) by spontaneous activation by background calcium, it will rapidly, regeneratively rise until most RyRs are open. Once this happens, casr will fall, raising the CICR threshold until it rises above nopen. At this point, CICR is no longer regenerative, and nopen will fall because of independent, spontaneous closure of RyRs. This will allow the depleted JSR to begin refilling, so the threshold will fall again, following nopen downward. If the falling threshold remains above the live value of nopen, nopen will converge to the stable, closed fixed point and remain there until a spontaneous opening triggers another spark.
The condition just described corresponds to the green curves in Fig. 7, which shows the coupled dynamics of the deterministic model Eqs. 1.5 and 1.6 (nopen in A, casr in B, and the phase–plane plot of nopen vs. casr in C). Those curves were computed with Dup = 0.1, a value that corresponds to the restitution time of 90 ms observed by Ramay et al. (2011). Suppose, however, we assume Dup = 0.7, the value derived from modeling the FSR as a fine, random network. In this case, our simulations revealed a rather complex behavior illustrated by the red curves in Fig. 7. What has happened is that the more rapid refilling of the JSR allows the falling CICR threshold to catch up with the falling nopen, so that CICR again becomes regenerative. The system can “latch up,” becoming permanently trapped in state with a small number of RyRs open, with their release flux balanced by steady refilling of the JSR, and the resulting cleft calcium sufficient to maintain their (average) open state. As shown in the phase–plane plot (Fig. 7 C), this indicates the creation of a new, stable fixed point of the combined system, onto which the trajectory converges instead of returning to the fully closed state. This corresponds, in the toy model, to the “metastable spark” observed in the stochastic model. As discussed further in the Appendix (especially Fig. A1), the presence of this new fixed point is a steady-state phenomenon, dependent on the relationship between the steady-state values of nopen and the steady-state depletion of the JSR, which depends on Dup but not on Casq. However, whether a metastable spark actually forms depends critically on the relative rates of the various processes. Relatively small changes in parameters can affect whether the trajectory is captured by the new equilibrium point.
In the more realistic, stochastic model, this trapped state is not permanently stable. Random fluctuations allow the system point to “leak” out of the attraction basin of the equilibrium point, permitting metastable sparks to extinguish slowly by a process of stochastic attrition, giving the long tail on the histogram of spark duration. Fig. 7 D shows phase plots of the ensemble average nopen and [Ca2+]jsr of 10,000 stochastic sparks in stable and metastable regimes. In the metastable regimen, sparks pause near a “latch point” before slowly returning to the resting state by stochastic attrition. This is a fundamental difference between deterministic and stochastic (more realistic) behavior.
What the toy model shows is that spark metastability is a kinetic matter, depending on the quantitative relationship among the rate of refilling of the JSR and the gating rates of the RyR. Anything that makes the JSR refill more quickly, the RyR more sensitive to calcium (including removal of regulation by lumenal calcium), or the RyR closure rate slower (open dwell time longer) will favor the metastable spark regimen. This is a matter of practical concern, because the prolonged release from metastable sparks may permit CICR to propagate to adjacent couplons (not considered in this paper) causing diastolic calcium release events that may result in triggered arrhythmias. In fact, mutations of Ryr2 and CASQ2 having these kinetic effects are known to cause catecholaminergic polymorphic ventricular tachycardia in humans and animal models (Jiang et al., 2004; Cerrone et al., 2005, 2007; Liu et al., 2006; Chopra et al., 2007; Paavola et al., 2007). Likewise, the arrhythmias that may accompany heart failure, and ischemic and reperfusion injury, may be explained by changes (phosphorylation, oxidation, nitrosylation) that affect gating of the RyR in ways leading to spark metastability (Eisner et al., 2009). The conditions for spark metastability are also critically dependent on the diffusion coefficient of free calcium in the dyadic cleft, as discussed below. This variable is not generally amenable to modification, but its value in vivo is poorly known.
“Treating” spark metastability
Treatments that redress the kinetic balance should be helpful for such arrhythmias. Fig. 8 shows spark extinction time histograms (log ordinate scale) for a model with RyR gating rates taken from a two-state approximation to the lipid bilayer data of Guo et al. (2012) (see supplemental text for a discussion of the issues in making this estimate) and assuming Dup = 0.15. As seen from the black curve, these parameters (taken from observations in different kinds of experiments) actually give sparks that are somewhat metastable. The green curve shows the effect of simultaneously doubling the RyR opening and closing rate constants, keeping the po as a function of [Ca2+]cleft unchanged. This kinetic change substantially reduces the tail of metastable sparks. The red curve shows the additional effect of doubling the concentration of Casq, thereby slowing JSR depletion and refilling. Increasing Casq delays the peak of the histogram as expected, but the tail actually moves leftward as the refilling rate decreases compared with the RyR closure rate, reducing spark metastability. This offers a possible explanation for the observation that drugs, such as carvedilol and flecainide, that increase the RyR closure rate (decrease mean open time) have been effective in preventing calcium waves and catecholaminergic polymorphic ventricular tachycardia in animal models (Watanabe et al., 2009; Hilliard et al., 2010; Zhou et al., 2011).
The paradox of restitution
As demonstrated by Picht et al. (2011) and confirmed by our simple spherically symmetric spark/blink computations (see above and supplemental text), the observed rates of blink recovery (∼200 ms in rabbit ventricle) and spark restitution (∼90 ms in rat using repeating sparks generated by micro-dose ryanodine; Sobie et al., 2005; Ramay et al., 2011) are not easily accounted for by simple models of calcium diffusion in the FSR lumen. Possible explanations are the presence of tortuosities or constrictions in the connection of the JSR to the FSR, or an error in the estimation of the intra-lumenal calcium diffusion rate. Picht et al. (2011), however, interpreted this discrepancy as evidence that there is a continuing tail of calcium release during the recovery from a blink, and that the time constant of recovery is limited by the decay of this release rather than by the rate of refilling from the FSR. Our analysis of spark metastability above makes this explanation problematic. If [Ca2+]SR recovers while some RyRs are still open, CICR will not extinguish: the spark would be metastable. On the assumption that metastable sparks are not a normal condition, we take this as evidence against Picht’s interpretation and favor the assumption that greater limitations on intra-SR diffusion are actually present. However, in the interest of open-mindedness, we must at least consider the hypothesis that “normal” spark/blink pairs are, in fact, metastable. For sparks that are terminated by a JSR depletion signal, the time-to-peak of spark fluorescence depends almost entirely on JSR volume, the number of RyRs recruited, and the resistance to efflux from the cleft. A metastable spark would only be recognized by a trailing “ember” (González et al., 2000) because of continuing release through open RyRs, fed by calcium returning from the FSR. Fig. 9 shows an overlay of 50 metastable sparks (F/F0 at the spark center). For these simulations, refilling rate was low but metastability was produced by raising the calcium sensitivity of the RyR. The resulting embers, although clearly visible in the noise-free simulation, might be difficult to detect under the conditions of confocal recording. Note also that the embers tend to be oscillatory. So called “sub-sparks” have in fact been noted after major sparks (Santiago et al., 2013) but have usually been attributed to recruitment of sub-clusters of RyRs. In our simulations, couplons are single crystalline arrays of RyRs, without satellite groupings, so oscillatory embers are a purely dynamical phenomenon (see Videos 1 and 2).
Spark amplitude distribution
In chemically skinned cat atrial cells, Shkryl et al. (2012) measured the distribution of amplitudes and rise times for a population of in-focus sparks, using a rapid two-dimensional confocal scanner combined with piezo z movement to locate the center of the spark in the z direction. Their data (replotted in Fig. 10 A) show that amplitudes are broadly distributed, with a long tail in the direction of higher amplitudes. For comparison, Fig. 10 B shows simulated spark amplitude distributions for several couplon sizes. For these simulations, we made the assumption that JSR volume is proportional to the number of RyRs in the couplon. That is, we assumed that the cell contains 1.4 × 106 RyRs (Bers and Stiffel, 1993) and that the whole-cell JSR occupies 0.35% of cell volume, partitioned equally among couplons. The larger couplons therefore produce larger sparks mainly because they have access to larger stores of JSR calcium. However, for any given couplon size, the amplitude distribution is rather narrow, with a sharp upper limit constrained by JSR contents. Fig. 10 C shows that for fixed couplon size, the amplitude distribution also depends on the RyR opening rate constant ko, which determines how many and how rapidly RyRs will be recruited and, therefore, the speed with which the JSR contents will be released. Even if RyR sensitivity is reduced to the point that sparks become highly stochastic, the amplitude distribution preserves the sharp upper cutoff. We conclude that the long tail on the spark amplitude distribution observed by Shkryl et al. (2012) cannot be accounted for by a homogeneous distribution of couplons with spark termination caused by local depletion. To account for the data would require either heterogeneity of couplon sizes or the presence of more complicated events involving local propagation of CICR between couplons. This is further confirmed by the joint distribution of amplitudes and rise times (Fig. 11). The observations show a positive correlation between rise time and amplitude. For simulated sparks, the relationship is negative for the fundamental reason that a longer rise time implies a slower release of approximately the same amount of JSR calcium.
Residual JSR calcium
Fig. 12 A shows the distribution of the nadir of JSR free calcium during a spark, with or without lumenal modulation of RyR gating. The fraction of calcium remaining in the JSR after a spark is quite low, as has been found in other models (Sobie et al., 2002; Ramay et al., 2011; Laver et al., 2013). This is at variance with blink observations showing only ∼50% depletion, as well as with the common idea that JSR calcium falls to a “threshold” of ∼50% at which point lumenal deactivation terminates the spark.
Although the limited resolution of confocal microscopy probably explains the relatively shallow blinks (Kong et al., 2013), we attempted to determine whether different conditions might allow spark termination with a higher true JSR calcium nadir. Reducing the RyR open time to 1 ms to allow RyRs to close rapidly once the CICR threshold is reached, and reducing the RyR calcium sensitivity, can increase the average residual calcium. However, the nadir of JSR calcium did not approach 50% unless RyR sensitivity was so low that sparks became very stochastic, with low amplitude and a very broad multimodal distribution of residual JSR calcium rather than a “threshold” (Fig. 12 B). In this condition, only 7% of spark attempts ignited a detectable spark. The difficulty in reproducing the hypothesized lumenal threshold behavior is caused by the fact that a substantial number of open RyRs will deplete the JSR very rapidly on a time scale comparable to RyR gating. This is a fundamental constraint caused by the assumed ratio of RyRs to JSR volume and the magnitude of the RyR unitary current. We varied the parameters of cytosolic and lumenal calcium regulation of RyR gating in an attempt to determine whether there was a regimen in which lumenal “threshold” behavior was observed. The only manipulation that was able to recreate such a threshold was to make the RyR closure rate (i.e., open dwell time) strongly and very nonlinearly dependent on lumenal calcium, effectively “slamming shut” the RyRs when lumenal calcium falls to the threshold (Fig. 12 C). As there is no experimental evidence for such strong lumenal regulation of RyR closure rate, our simulations support the belief that the apparent 50% residual calcium in “blinks” is caused by limitations of confocal imaging of the JSR terminal, which is at the resolution limit of light microscopy. The issue remains unresolved, however, as there is strong experimental evidence (namely, agreement between depletion measured in blinks and in more global images) suggesting that the limiting depletion is well measured locally (Zima et al., 2008).
Extent of RyR recruitment
Somewhat surprisingly, for the larger couplon sizes not all RyRs open even during a robust spark. Fig. 13 A shows histograms of the number of activated RyRs for several RyR opening rate constants, as in Fig. 12 C. The explanation is that lumenal calcium falls at a rate comparable to the rate at which activation of RyRs propagates across the couplon (Fig. 13 B), so that CICR regeneration may fail while many RyRs have not yet been opened. The propagation of RyR activation is shown in Videos 1 and 2.
Diffusion of free calcium in the dyadic cleft
The resistance to exit of calcium from the dyadic cleft depends on the local diffusion coefficient of free calcium (DCa), which is not directly known. In free solution, DCa depends on temperature, having a value of 0.748 µ2/ms at 23°C and 1.05 at 37°C (Hrabetová et al., 2009). This value is reduced by viscosity and “tortuosity” in cytosol; values of 0.25–0.35 have been assumed in spark modeling studies (Smith et al., 1998; Soeller and Cannell, 2002; Laver et al., 2013). The dyadic cleft is largely occupied by RyR protein, so it is likely that further tortuosity effects exist. Tadross et al. (2012) found that the calcium nano-domain in the vicinity of an L-type channel pore is 10 times stronger than predicted, implying a highly restricted diffusion environment (i.e., DCa ≈ 0.07). For most of our simulations, we used DCa = 0.15.
Fig. 14 shows the distributions of spark extinction times assuming RyR kinetics derived from the results of Guo et al. (2012) as described above, and using values of DCa of 0.7, 0.35, 0.15, or 0.07 µ2/ms. It is apparent that the degree of spark metastability is markedly increased by restricted calcium diffusion in the cleft; for the value of DCa implied by Tadross et al. (2012), 75% of sparks fail to terminate within 1 s. This problem cannot be dealt with by rescaling the model: cleft geometry is constrained by ultrastructure, JSR contents are constrained by ultrastructure and biochemical measurements of CSQ, and RyR kinetics and conductance are constrained by lipid bilayer studies, so CICR loop gain is determined once DCa has been specified. This places a premium on understanding the extent to which diffusion and channel properties are modified by the in situ environment of the cleft.
Calcium sparks in cardiac muscle are generally understood to be caused by locally regenerative CICR within a cluster of RyRs. How such regeneration, once ignited, is terminated is a matter of uncertainty and controversy. A commonly espoused idea is that local depletion of calcium in the JSR causes deactivation of RyRs–by a luminal calcium–sensing site on either the RyR itself or Casq–once JSR calcium falls below a critical threshold (Lukyanenko et al., 2001; Zima et al., 2008). We performed numerical experiments on a simplified model of stochastic CICR to test whether this mechanism can terminate sparks in the absence of any inactivation of RyRs. Using a minimal two-state RyR gating scheme, the initiation of sparks by a single RyR opening and their termination in response to JSR depletion could be reproduced. The model was deliberately simplistic to study the mechanisms in a reductionist manner. Therefore, we made no attempt to fit quantitatively the experimental data. To do so would be unwise, in any case, because crucial parameters of RyR gating, diffusion within the SR, and buffering in the cytosol show wide discrepancies among studies.
The most important finding in our study is that the ability of depletion to terminate sparks is not absolute but depends on the kinetics of RyR gating and JSR refilling. When RyR calcium sensitivity is high, RyR open dwell time is long, and/or intra-lumenal diffusion of calcium in the SR is rapid or buffering in the JSR is weak, sparks can fail to terminate cleanly. Instead, they enter a period of prolonged, stuttering, often oscillatory, low grade release. Such release is ultimately terminated randomly by stochastic attrition, leading to a long tail on the distribution of spark durations. We refer to this phenomenon as “spark metastability.”
Such prolonged RyR opening has the potential to trigger calcium release from adjacent couplons, resulting in propagated diastolic CICR or waves, which can trigger life-threatening arrhythmias. Spark metastability might be difficult to recognize when it is produced by high RyR sensitivity with normal refilling of the JSR, as the trailing calcium release from the depleted JSR might be inconspicuous. Zima et al. (2008) demonstrated extremely prolonged sparks when RyRs were partially blocked by tetracaine. They interpreted this as evidence that refilling of the JSR was rapid enough to prevent [Ca2+]JSR from being depleted to the “threshold” at which release would be terminated. The observation is similar to our spark metastability, except that release flux appeared to remain constant for hundreds of milliseconds, as though only a single channel were open or RyRs gated synchronously. This phenomenon appears to be specific to tetracaine, as no similar prolonged openings were seen when RyR unitary current was reduced by Tris (Guo et al., 2012), which should have similar effects if the tetracaine effect is caused by a reduced efflux of JSR calcium through partially blocked RyRs, as proposed by Zima et al. (2008). Long-lasting openings of the kind seen by Zima might well be the result of a more complicated gating scheme than the minimalist one we deliberately restricted ourselves to. It is possible that tetracaine produces long-lasting open states or subconductance states (as ryanodine does) that have not been studied, that there is some form of inactivation that appears in vivo although not seen in bilayer and, most intriguingly, that cooperativity caused by allosteric interactions might result in collective modes. These possibilities can be explored in future work using the same model structure but open up a vast territory of parameters and gating schemes that are not considered here.
Somewhat unexpectedly, the ability of JSR depletion to terminate sparks did not require the modulation of RyR gating by lumenal calcium. In the absence of such modulation sparks still terminated. This indicates that the interruption of positive CICR feedback by decreased RyR flux is a sufficient and probably dominant mechanism of spark termination. A similar conclusion has been reached, indirectly, from the effects of reducing RyR unitary current with Tris (Guo et al., 2012), and from recent modeling studies by Laver et al. (2013) using a model of a couplon with slightly more realistic geometry than ours. Even when lumenal modulation of RyR opening rate was present, we found little evidence for the existence of a deterministic “threshold” of SR calcium at which sparks would terminate. Unless RyR closure rate was made strongly dependent on lumenal calcium, sparks had a very low nadir of JSR calcium. Similarly strong depletion of JSR calcium has been found in other models (Sobie et al., 2002; Ramay et al., 2011; Laver et al., 2013). This suggests that the shallow depletion during “blinks” is an artifact of blurring by confocal microscopy. This has been demonstrated recently by Kong et al. (2013), who used a deterministic spark diffusion model (very similar to the one we used in Fig. S3) to fit observed spark/blink pairs, without assumptions about the underlying release mechanism. They found that the optimal fit required marked JSR depletion (masked by blurring), early decay of release flux, and RyR opening lasting well past the peak of fluorescence, as produced in our stochastic simulations.
The role of JSR depletion in terminating sparks, and the associated phenomenon of spark metastability, places a premium on understanding the rate of refilling of the JSR. Diffusion in the SR lumen has been measured by Picht et al. (2011) using fluorescence photobleaching recovery, obtaining a value of 0.06 µm2/ms, which is an order of magnitude less than free calcium diffusion in aqueous solution, and approximately one fifth of the values usually assumed to pertain in cytosol. Their model, which assumed that there is a concentration of 27 mM of calcium-binding sites on Casq in the JSR, estimated JSR refilling rates at least three times faster than observed blink recovery rates. Our naive computation, assuming a random FSR network extending directly to the JSR, gave refilling rates several times faster yet. Picht et al. (2011) rationalized this by assuming that there is continuing release during the blink recovery phase. As we have shown, this can only be true if sparks are metastable. Alternatively, one must consider the presence of an additional diffusion resistance between the FSR and JSR, or a much higher concentration of buffer in the JSR. Brochet et al. (2005) estimated in rabbit ventricular myocytes that there is an average of 4.5 discrete, narrow connections attaching the JSR release terminal to the FSR network. The concentration of Casq in the JSR is also poorly known, with estimates ranging from 10 to 120 mM of Ca-binding sites (Murphy et al., 2011).We therefore assumed the presence of a phenomenological diffusion resistance connecting the JSR with the distant FSR pool, characterized by a refilling rate parameter Dup. In the simplest form of the model, we treated this as a single resistance between the JSR and the infinite pool of SR calcium far from the spark. Because Picht et al. (2011) showed in their model that SR depletion during a blink can extend for several microns, we also constructed a version of the model in which the local FSR is represented as a long tube (Fig.1 B) that is able to supply some of its calcium to the JSR during the course of a single spark. This resulted in somewhat larger spark amplitudes as expected (not depicted) but did not change the qualitative behavior of the model.
As discussed in the supplemental text, we attempted to estimate the parameters of our two-state RyR gating scheme from lipid bilayer measurements made in “cell-like salt solutions” (Guo et al., 2012). These resulted in sparks that were somewhat metastable, so most of our simulations were done using lower values of the opening rate constant ko (equivalently, higher EC50 of the open probability curve). Laver et al. (2013) made stochastic simulations of sparks based on measured parameters and obtained sparks that terminated promptly. This results from two factors. For unknown reasons, their lipid bilayer data in rat showed an EC50 an order of magnitude higher than that found by Guo et al. (2012). Second, although most of the geometry of their sarcomere was based on observations, they assumed that the JSR filled from the FSR through a single 30-nm “nexus,” which has an effect similar to our imposition of an empirically constrained low value of Dup.
By studying the analogue of spark metastability in a simple, deterministic “toy” model of the couplon, we established that spark metastability is fundamentally a kinetic phenomenon, determined by the relationship between the refilling rate of the JSR, which controls spark refractoriness, and the gating rates of the RyR, which determine how quickly the open probability of the RyR diminishes in response to JSR depletion, preventing re-ignition. Because uncontrolled CICR can have life-threatening consequences, we considered whether it might be possible to “treat” spark metastability by altering this kinetic balance. To demonstrate the kinetic nature of the effect, we increased both the opening and closing rate constants of the RyR by the same factor, leaving the po as a function of calcium unchanged. This greatly reduced spark metastability (Fig. 8), which was further improved by slowing JSR refilling by increasing Casq from 30 to 60 mM. Zhou et al. (2011) studied 14 β blockers used to treat ventricular arrhythmias in heart failure and found that only one, carvedilol, could suppress calcium waves. They showed that it acted by decreasing the open dwell time of the RyR, which is equivalent to increasing the closing rate kom. It had no effect on the opening rate constant (which might normally be thought of as the “calcium sensitivity” of the RyR) and actually increased the total number of openings. However, because the closing rate controls the rate at which RyR openings respond to JSR depletion, as well as the rate of stochastic attrition, increasing the closing rate has the most effect on spark metastability, explaining why carvedilol is most effective in treating calcium-overload arrhythmias (Poole-Wilson et al., 2003). Flecainide, another antiarrhythmic drug that reduces RyR open time, has also been shown to be uniquely protective against calcium waves and triggered arrhythmias (Watanabe et al., 2009; Hilliard et al., 2010).
An analysis similar to our toy model was presented by Hinch (2004), who examined the asymptotic properties of a simplified Markov model of sparks, leading to quasi-deterministic analytic results parallel to what we have done in the “toy.” Hinch examined a somewhat different regimen (basically that of Sobie’s “sticky cluster” model; Sobie et al., 2002), which favors bistable behavior (flat top sparks) with a narrow stochastic transition between open and closed. The phenomenon that we call “spark metastability” was demonstrated, using a null-cline analysis similar to our Fig. A1 (Appendix), specifically as an effect of removing allosteric coupling, although the generality of the effect as a property of relative kinetic rates was not explored.
A model is actually most informative when it is unable to fit the data. Figs. 10 and 11 show that in at least one model system–chemically skinned cat atrial myocytes–the distribution of in-focus spark amplitudes has features that cannot be explained by our model. The fact that the observed amplitude distribution has a long tail that is positively correlated with spark rise time cannot be accounted for by a model in which spark termination is caused mainly by JSR depletion in a homogeneous set of couplons. One can speculate that this might be evidence of heterogeneity of couplon properties, or, more radically, that events counted as “sparks” actually involve local propagation to more than one couplon. This latter possibility seems especially plausible because the atrium, lacking t-tubules, depends on radial propagation of CICR to excite the center of the fiber. In service of this propagation, there may be small, nonjunctional RyR clusters with spacings as close as 100–150 nm (Franzini-Armstrong, C., personal communication) that might appear as a single “spark” source. This invocation of heterogeneity may also be appropriate to account for another observation by Shkryl et al. (2012) that is not readily consistent with any of the termination mechanisms discussed here, namely, that sparks of different rise time are essentially superimposable (in fluorescence and its rate of rise) up to their peak time.
Summary and future work
By means of numerical simulations, we have demonstrated that it is possible to achieve termination of calcium sparks in a model with RyRs having only two states, with no inactivation or adaptation process. Termination results from local depletion of JSR calcium, which is an inevitable consequence of the high conductance of the RyR channel and the architecture of the SR. Depletion terminates sparks by interrupting the positive feedback of CICR by reducing calcium flux into the dyadic cleft, and also by reducing the open probability of the RyR by removal of activation by lumenal calcium, although the latter mechanism is not essential. However, termination of sparks by these mechanisms is not absolute and can be greatly prolonged if the kinetic relationship among gating rates of the RyR and refilling of the JSR is unfavorable, a phenomenon we dub “spark metastability.” This has the potential to contribute to the development of life-threatening triggered arrhythmias. Spark metastability is critically dependent on the calcium-dependent opening and closing rates of the RyR, the resistance to diffusion within the SR lumen, the concentration of Casq in the JSR, and the diffusion coefficient of free calcium in the dyadic cleft, none of which is known with precision. In particular, more evidence is needed regarding restriction of diffusion within the dyadic cleft, which our simulations (Fig. 14) show would have a profound effect.
Using the best estimates of these factors, sparks from couplons with ∼50 RyRs may, in fact, be metastable. This indicates either that metastable sparks do occur in vivo, or that further experimental work is needed to better determine the conditions at the couplon. Particular uncertainty attaches to the rate of refilling of the JSR by diffusion within the SR lumen, as observed rates are severalfold slower than might be expected from simple models. Restrictions to diffusion of calcium in the cleft need further study, as recent estimates are incompatible with spark termination in large couplons.
The simulations shown here assumed a cooperativity of 2 for calcium activation of the RyR, as that is usually found in lipid bilayer in “cell-like” solutions with high lumenal-side [Ca2+]. However, it is difficult to explain the low resting spark rate under these conditions. Our qualitative findings are not altered by using higher cooperativity for RyR opening (not depicted), but the precise conditions for spark metastability may be different. One possible mechanism to explain high cooperativity in vivo would be the existence of allosteric interactions among RyRs in dense clusters (Stern et al., 1999). We have not included allosteric interactions in spark simulations because evidence for its existence is sparse, but it needs to be studied in future simulations, together with the fact that couplons have been found to consist of smaller sub-clusters (Baddeley et al., 2009).
Finally, we must mention the most important idealization in our model: the fact that sparks are treated as isolated and independent. It is conventional to say that sparks are the elementary events of excitation–contraction coupling, from which whole-cell calcium transients are formed by superposition. In fact, however, calcium release from one couplon will influence its neighbors in at least three ways: (1) cytosolic calcium will diffuse to activate, or at least bias, RyRs; (2) SR calcium will be drawn away from other couplons by diffusion within the lumen; and (3) SR calcium will be increased by diffusion through the cytosol followed by uptake by SERCA. Each of these effects has different spatio-temporal characteristics. The concept of CICR metastability needs to be extended to the whole cell, which is a hierarchical array of RyRs, a kind of couplon of couplons. Such extension goes beyond the scope of this paper.
Construction of the “toy” model
We first estimate the average free calcium in the dyadic cleft when a calcium current i is entering the cleft. We assume that the distribution of calcium in the cleft is in instantaneous steady state. The cleft is treated as a cylinder of radius R and (small) height h and analyzed in two dimensions. At the edge of the cleft, a boundary value of ca0, equal to bulk cytosolic [Ca2+], is assumed. This is a crude approximation but good enough for our qualitative purposes. The calcium current i entering the cleft via RyRs is taken to all be deposited at the center.
In this approximation, free calcium as a function of radius will satisfy the two-dimensional Laplace equation in polar coordinates, with circular symmetry, except at the origin, which is singular because of the delta-function source:(A1.7)
The solution satisfying the boundary condition is easily verified to be:(A1.8)
Differentiating and setting r = R, we find that the flux of calcium per unit area exiting at the edge of the cleft is given by:(A1.9)from which the constant a can be determined by requiring that the total calcium flux out of the cleft be equal to the calcium entering via the current i. This gives the solution:(A1.10)where F is the faraday.
Integrating the above expression over the area of the cleft gives the average [Ca2+] in the cleft, which we will take as the calcium seen by RyRs on their cytosolic face:(A1.11)
Note that this expression depends on the height of the cleft but not on its radius. If we had assumed that the source current were uniformly distributed over the area of the cleft, we would have obtained a parabolic expression for cleft calcium, giving an average value differing from the above by a factor of two, which is immaterial for our qualitative purposes. Recall that, in the full, stochastic model, the exact distribution of cleft calcium is calculated dynamically at the position of each RyR.
The current i is taken to enter through nopen RyRs, each of which has linear permeation proportional to the free calcium gradient across the pore, such that the unitary current under resting conditions (casr = casr0, <ca> = ca0) is iu0:(A1.12)
Eliminating <ca> between these two equations gives:(A1.13)
The total calcium concentration in the JSR, including instantaneous buffering by Casq, is:(A1.14)
casrt changes caused by calcium leaving through RyRs and being replenished from the distant SR pool through a single diffusion resistance:(A1.15)
Performing the differentiation and solving for the rate of change of free calcium in the JSR gives:(A1.16)
Inserting the above-derived expression for i gives the casr dynamical Eq. 1.5 shown in the main text:(A1.17)
Setting the left-hand side to zero and solving for casr gives the steady-state value that JSR calcium would achieve if nopen RyRs were kept open indefinitely, which is plotted in Fig. 2, for various values of Dup:(A1.18)
We assume that the continuous variable nopen changes in the same way as the ensemble average for RyRs activated by a cytosolic-side calcium <ca>:(A1.19)where hry is a Hill coefficient giving the cooperativity of RyR activation (taken to be 2 in the numerical examples in the text).
Inserting the above-derived expression for <ca> gives the dynamical Eq. 1.6 for nopen shown in the main text:(A1.20)
Phase–plane stability analysis of the toy model
Here, we further clarify the steady-state and kinetic aspects of spark metastability using the deterministic toy model. The model provides an approximation to spark metastability as trapping of the system trajectory on a stable fixed point, an aspect of the steady-state behavior of the system. Fig. A1 shows how to locate such a point by a graphical analysis of the phase plane of nopen versus casr for stable (A) and metastable (B) parameter values (as in Fig. 7, but plotted with logarithmic ordinate for clarity). In each panel, the dashed thin black line shows the curve (nullcline) along which dnopen/dt is zero. For points above the curve, nopen will increase, whereas for points below, it will decrease. This curve is the same in both panels; its location does not depend either on Dup or on the absolute scale of RyR gating rates or the presence of Casq. The thin magenta line shows the nullcline along which dcasr/dt vanishes; i.e., it plots the steady-state depletion of the JSR. This is the same information as in Fig. 2 but plotted as a function of nopen for fixed Dup. The location of this curve depends on Dup but not on RyR gating rates or Casq. A fixed or stationary point of the full system will exist wherever these two curves intersect. For the parameters in Fig. A1 A, where Dup = 0.1, there is no intersection, so no “latch point” exists and no spark metastability is possible. The system trajectory (Fig. A1 A, green curve; same as in Fig. 7 C except for the logarithmic axis) returns exponentially to the resting (nopen = 0) state. In Fig. A1 B, Dup = 0.7. This raises the magenta curve so that it now intersects the dashed nopen stationary curve in two places, creating a pair of fixed points, stable (bottom right) and unstable (top left). The de novo appearance of these equilibria as a function of a parameter is an example of a saddle bifurcation. However, the existence of a stable fixed point does not guarantee spark metastability, which is determined by whether the trajectory is actually captured by the point. That is a matter of kinetics. For the particular parameters in Fig. A1 B, the red trajectory is captured and spirals into the fixed point where it will remain. However, increasing Casq from 27 to 120 mM, which does not affect static stability (i.e., the location of the equilibrium point), allows the trajectory (blue curve) to circumvent the equilibrium point and return to the resting state. The same effect can be produced by increasing both the opening and closing rates of the RyR proportionately (not depicted). The possible role of Casq buffering to stabilize sparks was pointed out by MacLennan and Chen (2009), but they interpreted it in terms of steady-state buffering preventing SR calcium from reaching a hypothetical threshold for store-overload release.
We would like to thank Professor Thomas Shannon of Rush University for sharing data and custom software, Professor Mark Cannell of University of Bristol, UK, for providing technical assistance with the FACSIMILE software, and Professor Michael Fill of Rush University for helpful discussions.
This work was supported by the Intramural Research Program of the National Institute on Aging, National Institutes of Health (NIH).
This study used the high performance computational capabilities of the Biowulf Linux cluster at the NIH.
Richard L. Moss served as editor.
- Abbreviations used in this paper:
- free SR
- junctional SR
- Submitted: 22 May 2013
- Accepted: 1 August 2013
This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 3.0 Unported license, as described at http://creativecommons.org/licenses/by-nc-sa/3.0/).