The role of cardiac ryanodine receptor (RyR) gating in the initiation and propagation of calcium waves was investigated using a mathematical model comprising a stochastic description of RyR gating and a deterministic description of calcium diffusion and sequestration. We used a one-dimensional array of equidistantly spaced RyR clusters, representing the confocal scanning line, to simulate the formation of calcium sparks. Our model provided an excellent description of the calcium dependence of the frequency of diastolic calcium sparks and of the increased tendency for the production of calcium waves after a decrease in cytosolic calcium buffering. We developed a hypothesis relating changes in the propensity to form calcium waves to changes of RyR gating and tested it by simulation. With a realistic RyR gating model, increased ability of RyR to be activated by Ca2+ strongly increased the propensity for generation of calcium waves at low (0.05–0.1-µM) calcium concentrations but only slightly at high (0.2–0.4-µM) calcium concentrations. Changes in RyR gating altered calcium wave formation by changing the calcium sensitivity of spontaneous calcium spark activation and/or the average number of open RyRs in spontaneous calcium sparks. Gating changes that did not affect RyR activation by Ca2+ had only a weak effect on the propensity to form calcium waves, even if they strongly increased calcium spark frequency. Calcium waves induced by modulating the properties of the RyR activation site could be suppressed by inhibiting the spontaneous opening of the RyR. These data can explain the increased tendency for production of calcium waves under conditions when RyR gating is altered in cardiac diseases.
In mammalian cardiac myocytes, calcium release from the SR through the RyRs is the source of calcium for contraction. RyRs are calcium release channels clustered at SR cisternae that are a part of the tubulo-reticular junction (Franzini-Armstrong et al., 1999). Their activity is dynamically regulated by cytosolic as well as luminal calcium concentration (Fill and Copello, 2002; Laver, 2007a; Györke and Terentyev, 2008; Tencerová et al., 2012), and modulated by Mg2+ ions (Rousseau et al., 1986; Meissner and Henderson, 1987; Ashley and Williams, 1990; Zahradníková et al., 2003, 2010; Laver, 2007b; Laver and Honen, 2008; Diaz-Sylvester et al., 2011; Tencerová et al., 2012).
During normal heart function, calcium release through RyRs is triggered by activity of voltage-dependent plasmalemmal calcium channels, a process named “Ca2+-induced Ca2+ release” (Fabiato, 1985; Bers, 2001), and occurs as a spatiotemporal summation of elementary calcium release events: calcium sparks (Cheng et al., 1993; Cannell et al., 1994). Under many pathological conditions, such as calcium overload (Cheng et al., 1993, 1996), heart failure (Yano et al., 2005), and genetic diseases that alter RyR function (Paavola et al., 2007) or calcium handling in the SR (Knollmann et al., 2006), spontaneous calcium release that may play a role in arrhythmias (Lakatta, 1992; Knollmann et al., 2006; Paavola et al., 2007) occurs in the form of calcium waves (Cheng et al., 1993; Lukyanenko and Gyorke, 1999) that are also a result of spatiotemporal summation of calcium sparks (Cheng et al., 1993, 1996). Although the occurrence of calcium waves in pathological situations has been ascribed to changes in RyR gating (Fernández-Velasco et al., 2009), the relationships between RyR activity and initiation and propagation of calcium waves are not clearly understood.
In this work, we propose a theoretical basis for understanding formation of calcium waves under diastolic conditions, based on a mechanistic model of the calcium spark (Zahradníková et al., 2010; Zahradníková and Zahradník, 2012). The theoretical predictions were verified using computer simulation of the stochastic process of calcium spark activation in a modified fire–diffuse–fire (FDF) model suggested here as an extension of the previous FDF models (Keizer et al., 1998; Coombes and Timofeeva, 2003). The results show that calcium waves can form and propagate at low frequency of spontaneous calcium sparks and in the physiological range of calcium release flux if the calcium dependence of spontaneous calcium spark frequency (which is a function of the RyR open probability) is sufficiently steep, or if the number of open RyRs in spontaneous calcium sparks is sufficiently large. Changes in RyR gating that increase these properties of spontaneous calcium sparks thus increase the propensity to form calcium waves.
MATERIALS AND METHODS
The aim of the presented theoretical approach was to quantitatively link the molecular properties of the RyR, namely its open probability governed by a gating mechanism, to the experimentally observed phenomena such as calcium sparks and waves occurring at the cellular level. Such a theoretical basis may provide a useful tool for understanding complex mechanisms of calcium signaling in cardiac myocytes and for interpretation of the experimental data at different levels, and it may also suggest several predictions of the system behavior under various conditions.
RyR activation and properties of spontaneous calcium sparks.
Open probability of the RyR is an experimentally accessible quantity, and it may be expressed as a mathematical function of various factors, particularly the concentrations of calcium and magnesium ions. We used the allosteric homotetrameric gating model (aHTG model) of RyR gating (Zahradníková et al., 2010) that was shown to quantitatively describe both the properties of single RyR channels and the frequency and amplitude distribution of spontaneous calcium sparks in cardiac myocytes (Zahradníková et al., 2010; Janíček et al., 2012; Zahradníková and Zahradník, 2012). The steady-state open probability PO can be expressed as:(1)where KCa and KMg are the Ca2+ and Mg2+ dissociation constant of the ion-free closed-state C00; fCa and fMg are the allosteric factors coupling Ca2+ and Mg2+ binding to channel opening; KO00 is the equilibrium constant of the O00–C00 transition (inversely proportional to the spontaneous opening tendency of the RyR); KOL, KCL, and KCI are the dissociation constants of the L-mode states OL and CL and the I-mode state I; KI is the Mg2+ dissociation constant of the divalent inhibitory site; and [Ca2+] and [Mg2+] are the free cytosolic concentrations of calcium and magnesium. Parameter values are summarized in Table 1.
The RyRs are organized in clusters that form calcium release units (CRUs), discrete calcium release sites distributed throughout the myocyte (Chen-Izu et al., 2006; Soeller et al., 2007; Baddeley et al., 2009; Zahradníková and Zahradník, 2012). The calcium spark is initiated by the first stochastic opening of an RyR in a CRU. The frequency of spontaneous calcium spark generation is therefore directly related to the steady-state open probability of a single RyR (Zahradníková et al., 2010):(2)where ϕspark is the frequency of spontaneous sparks per CRU, trd denotes the duration of release, trp is the period during which the CRU remains refractory after calcium release, τo represents the mean RyR open time, and nRyR is the median number of RyR channels in the CRU (Zahradníková and Zahradník, 2012).
The average period between spontaneous sparks at a single CRU can be expressed as the sum trd + trp + tit, where the idle time tit is a period of time when the CRU is neither refractory nor releasing (Zahradníková et al., 2010). The idle time is associated with the average closed time τC of a single RyR and with the median number of RyRs in the CRU, nRyR (Kunze et al., 1985; Zahradníková et al., 2010; Zahradníková and Zahradník, 2012):(3)Then, the steady-state probability of finding a CRU in the active state can be determined as(4)An analogous expression describes the steady-state probability of finding a CRU in the refractory state:(5)The probability of a CRU being neither active nor refractory is then(6)The probabilities of finding a CRU in a particular state (active–refractory–idle) expressed by Eqs. 4–6 are important for setting the correct initial states of the CRUs in the model.
When a CRU becomes active, it turns for a certain time period to a source of calcium ions of a defined strength. The amplitude of calcium release flux from each active CRU is directly proportional to the number of open RyRs in this CRU during a spark. Because the RyR opening is a stochastic process, calcium flux amplitudes are not uniform, but rather have quantal character, as observed experimentally (Wang et al., 2004; Janíček et al., 2012). We have shown recently (Zahradníková et al., 2010; Janíček et al., 2012; Zahradníková and Zahradník, 2012) that the probability that exactly nO RyR channels are open during a spark, PnO, can be determined using the binomial equation:(7)where pO is the probability that an additional RyR channel will open during a spark, as a result of recruitment by the first open RyR in the CRU. The probability pO can be calculated from the open probabilities of RyRs with a defined Mg occupancy and from the probability, pdis, that a Mg2+ ion dissociates from the activation site of the RyR during the spark (Zahradníková et al., 2010; Zahradníková and Zahradník, 2012). The procedure for pO calculation is given in the supplemental text (Eq. S1). Besides its use in Eq. 7, the probability pO allows for determination of the average number of RyRs that open during a spontaneous spark, which is a helpful quantity characterizing the spark amplitude:(8)The theoretical framework described above makes it possible to characterize and analyze in a quantitative manner the steady-state properties of spontaneous calcium sparks in the absence of wave formation using analytical expressions, as well as the dynamic properties of spontaneous and evoked calcium sparks using computer simulations. These expressions are used to derive other practical quantities and relations for further analysis, namely the sensitivity of the frequency of spontaneous calcium sparks (∂CaϕSpark) and of the average number of open RyRs in spontaneous calcium sparks to changes in cytosolic Ca2+ concentration (see Methods and Results). The theoretical predictions are compared with simulation results obtained using the spatiotemporal model (see the next subsection and Methods) of calcium sparks and waves in a cardiac myocyte, as described in the Results.
Dynamics of calcium release and diffusion.
Calcium ions released from a CRU diffuse in the cytosol and bring about an increase of Ca2+ concentration in the vicinity of nearby CRUs. As the open probability of RyRs strongly depends on [Ca2+], activation of one CRU may trigger calcium release from a neighboring CRU. In this way, calcium sparks may propagate in space and form a calcium wave. As described in the Introduction, transition of calcium sparks to calcium waves is an important phenomenon, associated with impaired calcium handling that underlies some cardiac arrhythmias. Quantitative description of the processes leading to formation of calcium waves requires analysis of the stochastic dynamics of CRU activation simultaneously with rendering calcium diffusion from multiple sources.
By restricting the calculation to a single dimension representing a one-dimensional regular array of CRUs, and by dividing the time domain into discrete time intervals Δt during which the stochastic events in question may occur, it is possible to apply the stochastic generalization (Coombes and Timofeeva, 2003) of the FDF approach (Keizer et al., 1998) to CRUs composed of clusters of RyRs described by a realistic gating model (Zahradníková and Zahradník, 2012). Within the FDF framework, it is assumed that changes of the CRU state occur only at integer multiples of the time step Δt. After activation, the CRU produces a rectangular pulse of Ca2+ ions of duration trd and amplitude JCa. After termination of release, the CRU becomes refractory for the period trp. Both trd and trp are restricted to integer multiples of Δt. The amplitude JCa represents the physical flux of Ca2+ (in mol m−2 s−1). Calcium is released uniformly to both directions. The flux JCa is given by the unitary flux jCa, corresponding to the calcium flux through a single open RyR, multiplied by the number of open RyRs in the CRU. Diffusion of calcium and its uptake by sarco/endoplasmic reticulum calcium ATPase (SERCA) is described by the following equation (Coombes et al., 2004):(9)
In Eq. 9, nCRU is the number of CRUs; mi is the number of calcium sparks occurring at the ith CRU; is the time at which the jth calcium spark of the ith CRU occurred; u denotes Ca2+ concentration at time t and location r; is the number of RyRs open during the jth calcium spark at the ith CRU; jCa is the unitary calcium flux, i.e., the amount of Ca2+ (in moles) that flows through a unit area per unit time if a single RyR is open; τd is the time constant of calcium removal from the cytosol by the SERCA pump; DCa stands for the effective diffusion coefficient of calcium ions that reflects Ca2+ buffering by endogenous and exogenous buffers; and the vectors ri represent the positions of the sources. Spatial and temporal discretization of the sources is expressed via the Dirac delta function δ and the Heaviside step function H, respectively. This means that the sources localized at ri are activated at the moments and terminated after the release duration trd. Under the assumption that release occurs on a regular temporal lattice, the solution of Eq. 9 can be constructed using the Green’s function and the fast-Fourier transform (FFT) technique, as described previously (Coombes et al., 2004).
The times necessary for solution of Eq. 9 are determined by a stochastic approach by considering state transition probabilities of the CRUs during each time interval Δt. The probability that a particular CRU is activated (i.e., that the calcium source is switched on) corresponds to the probability that the particular CRU changes its state from non-active (idle) to a state with at least one RyR channel open during the time Δt. The probability that in a CRU containing nRyR channels at least one opens during Δt, PCRU, is equal to the negation of the probability of the complementary event (no channel opens in the CRU), i.e., it can be expressed as:(10)where PRyR is the probability that a single RyR channel opens during the interval Δt (Inoue and Bridge, 2005; Poláková et al., 2008; Zahradníková et al., 2010). This probability can be found as the negation of probability that a single RyR remains closed for the time Δt, which is (Smith, 2002), where k+ represents the rate constant of the transition to the open state. Considering the relationship between k+, single-channel steady-state open probability PO(Ca), and the average open time τo, the probability that a single RyR opens within the interval Δt is(11)The introduced theory and the derived equations were used to construct the simulation model for analysis of calcium sparks and calcium waves in a cardiac myocyte.
Sensitivity of spontaneous spark properties to changes in cytosolic [Ca2+].
Our aim was to find the relationships between the properties of RyRs, the properties of spontaneous calcium sparks, and the propensity to form calcium waves. To this end, we characterized spontaneous calcium sparks by their average frequency (ϕspark) and average number of open RyRs Because formation of calcium waves is related to the ability of spontaneous calcium spark to ignite a triggered calcium spark at a neighboring release site by elevating the local cytosolic calcium concentration, we have introduced two additional empirical quantities: the sensitivity of the frequency of spontaneous sparks and the sensitivity of the average number of open RyRs in spontaneous calcium sparks to changes in cytosolic [Ca2+]. On one hand, the theoretical value of these four quantities can be unequivocally calculated from the model of RyR gating and spark generation used in this study; on the other hand, these quantities can be experimentally measured and serve to characterize spontaneous calcium sparks empirically, without a need for the knowledge of the underlying RyR gating model.
To analyze the effect of model parameters on the sensitivity of the frequency of spontaneous sparks (ϕspark) to changes in cytosolic calcium concentration, the quantity ∂Caϕspark = ∂ϕspark/∂[Ca2+] was introduced. The analytical expression for ∂Caϕspark was derived by differentiating Eq. 2, into which Eq. 1 has been substituted, with respect to [Ca2+]. The differentiation was performed analytically using Mathematica (version 9; Wolfram Research), and the complete formula is given in the supplemental text (Eq. S7).
Similarly, the quantity was introduced to quantify the sensitivity of the average number of open RyRs in spontaneous calcium sparks to the changes in cytosolic [Ca2+]. Because calculation of requires the knowledge of pO (Eq. S1), which is not defined by an analytical expression, was determined as a function of cytosolic [Ca2+] for different values of parameters of the aHTG model by numerical differentiation of Eq. 8 using OriginPro software (OriginLab).
Spatiotemporal model of calcium dynamics.
Spatiotemporal dynamics of calcium sparks and waves was simulated using the FDF approach (Keizer et al., 1998) with stochastic generalization (Coombes and Timofeeva, 2003) on a refined time grid with realistic composition of RyR clusters in the CRUs (Zahradníková and Zahradník, 2012) and a realistic model of RyR gating (Zahradníková et al., 2010), as described in the Theory section.
Initial conditions and model parameters.
The simulation domain was represented by a one-dimensional array of 51 equidistantly distributed calcium release sites corresponding to a typical confocal scanning line. Each CRU clustered nRyR = 36 RyRs, and the distance between neighboring CRUs was set to d = 2 µm. The domain environment represents the cytosol with a constant basal calcium concentration (typically set to the diastolic value of 0.05 µM; De Young and Scarpa, 1987; Smogorzewski et al., 1998). Three possible states of individual CRUs are considered: closed state (non-active, idle), when all RyRs of the CRU are closed and the CRU is ready to be immediately activated; the active state, corresponding to the calcium spark, when at least one RyR of the CRU is open (Zahradníková et al., 2010; Janíček et al., 2012), and the CRU is a source of Ca2+ with a flux JCa for a period trd as described below; and finally the refractory state, after the active state and corresponding to the refractory phase of the CRU after a calcium spark (Sobie et al., 2005), when all RyRs of the CRU are closed, but the CRU cannot be activated for a time period trp. The simulation was initialized by randomly assigning active, closed, and refractory state to individual CRUs with probabilities conforming to the values (ssPspark, ssPidle, ssPref) given by Eqs. 3–6. The mean RyR open time was set to τo = 15 ms (Tencerová et al., 2012; Gaburjakova and Gaburjakova, 2014), and the probability PO was calculated from Eq. 1 with parameters listed in Table 1 (Zahradníková et al., 2010; Zahradníková and Zahradník, 2012). At the start of the simulation, the calcium concentration in the whole simulated space was equal to the basal [Ca2+]. Afterward, the local dynamically calculated [Ca2+] at each CRU (see below) was used to calculate its PCRU (Eq. 10) during the current time step. The cytosolic concentration of Mg2+ ions was assumed to be constant at 1 mM. The initially refractory CRUs were assigned to dwell in their states for a period set up by a random drawing from a uniform distribution on the interval (0, trp). The CRU parameters with their standard values are summarized in Table 2.
CRUs in the active state are calcium sources with a flux JCa, and the released Ca2+ ions diffuse throughout the domain with an apparent diffusion coefficient DCa that accounts for calcium binding to endogenous and exogenous buffers in the cytosol. Reuptake of Ca2+ ions from the cytosol into the SR by SERCA was represented by the time constant τd. The flux JCa is given by the number of open RyRs (nO) in the active CRU and the unitary flux jCa through a single RyR (JCa = nO × jCa). The actual value of nO for each active CRU was set randomly according to the binomial probability distribution as explained in the Theory section (Eq. 7). The value of nO was generated using the “binornd” function from the Statistical Toolbox of MATLAB (version 8.3 [R2014a]; MathWorks) as follows:(12)where the probability pO has been explained in the Theory section, and the procedure for pO calculation is given in the supplemental text (Eq. S1). The local calcium concentration at the given CRU at the time of activation was used for calculating pO.
The calcium release flux through a single open RyR (the unitary flux jCa) was assigned by the following procedure. The spatiotemporal distribution during three-dimensional radial diffusion of Ca2+ from a point source of a strength corresponding to a 0.4-pA current (Mejía-Alvarez et al., 1999; Janíček et al., 2012) in the absence of Ca2+ uptake by SERCA was calculated analytically (Crank, 1975; Simon and Llinás, 1985; Izu et al., 2001) and compared with the analytical solution of the one-dimensional model (Coombes et al., 2004). The jCa value was adjusted so that [Ca2+] at r = 2 µm and t = 15 ms was equal for both models for the standard value of the diffusion coefficient (DCa = 30 µm2 s−1). This procedure yielded the value jCa = 0.481 µmol m−2 s−1.
The time constant of Ca2+ sequestration by SERCA at low cytosolic [Ca2+], τdmin, was taken from Schuhmeier and Melzer (2004), and its calcium dependence was simulated according to Swietach et al. (2010) as described in the supplemental text. The resulting calcium dependence of τd is shown in Fig. S1. Importantly, the value of τd is only very weakly cytosolic calcium dependent at [Ca2+] ≤ 0.25 µM, which justifies the use of Eq. 9 for calculation of the spatiotemporal evolution of calcium concentrations.
Model parameters that relate to Ca2+ diffusion and sequestration are summarized in Table 3.
Simulation of calcium dynamics.
Determination of calcium concentration in space and time involves solution of Eq. 9 with the proper values of times for the on/off switching of individual CRUs. Although the solution of Eq. 9 is calculated analytically, the simulation has to be performed in a step-by step process, analogous to numerical temporal integration, caused by the stochastic nature of the CRU-state transitions. Hence, the simulation time was divided to equal time intervals Δt = 3 ms. To assure a steady state for the whole analyzed time period, the simulation was started at a time t = −trp with a randomly assigned initial state of the array of CRUs as described above. At every time interval (ti, ti + Δt) of the simulation, a uniform random number (interval from 0, 1) was generated for each of the CRUs that were not refractory and compared with the probability PCRU, calculated from the local value of cytosolic [Ca2+] at the given CRU at t = ti (Eqs. 10 and 11). If the random number did not exceed the value of PCRU, the CRU became a Ca2+ source with the flux nO × jCa for the duration trd and subsequently became refractory for the duration trp. The value of nO was generated randomly after the binomial distribution (Eq. 12). The spatial distribution of cytosolic [Ca2+] in the domain at t = ti was evaluated as the solution of Eq. 9, with the actually switched-on sources using the corresponding Green’s function and FFT in relation to the preceding solution at t = ti − Δt. The calcium sparks were analyzed in the interval from t = 0 to t = 10 s to yield the values of the observed frequency of spontaneous and triggered calcium sparks, and the observed average number of open RyRs per CRU in spontaneous and triggered sparks. Under control conditions, the maximum increase of cytosolic [Ca2+] produced at a neighboring release site by calcium release flux of 0.98 pA, corresponding to open RyRs (Fig. S2), was Δ[Ca2+] = 0.21 µM, and the increase exceeded 0.15 µM in the interval 18–30 ms after spark activation.
Analysis of calcium waves.
An analyzed calcium release event consisted either of a spontaneous solitary calcium spark (spark event), or of an initial spontaneous spark followed by one or more triggered sparks (i.e., wave event). A spark was considered a triggered spark if it was activated within an interval of 5–50 ms relative to the start of another spark at a neighboring CRU (corresponding to wave velocity of 40–400 µm s−1). Sparks occurring at an interval Δt < 5 ms or Δt > 50 ms were supposed to be independent spontaneous sparks. Under control conditions at Δt = 5 ms, the cytosolic calcium at a neighboring calcium release site is only negligibly increased (Fig. S2). Indeed, the number of sparks occurring at Δt < 5 ms was low and proportional to spontaneous spark frequency, justifying the first assumption. The second assumption also seems justified, considering that calcium waves in cardiac myocytes exhibit a wave velocity of ∼80 µm s−1 (Swietach et al., 2010), and that the cytosolic calcium elevation at a distance of 2 µm at 50 ms after spark activation has already declined to <25% of its peak value (Fig. S2 B).
The waves were characterized by their length, which was expressed as the average number of sparks per wave event, Termination of calcium waves was classified as caused by refractoriness of the release site 2 µm ahead of the wave front, as a result of collision with another calcium wave or because of propagation failure. Further details are presented in the supplemental text and illustrated in Fig. S3.
Software and data analysis
All derivations were obtained in Mathematica. Simulation model of calcium dynamics was implemented in MATLAB. All random numbers were generated using the Mersenne Twister algorithm (Matsumoto and Nishimura, 1998). Wave analyses were performed in MATLAB. Multiple linear fitting of the calcium dependence of spark and wave properties was performed in Mathematica. All other analyses were performed using the program OriginPro. Three simulations lasting 10 s were used to calculate the average values of the analyzed variables. Data are given as mean ± SE; fitted parameters are given with their standard error of the fit. Statistical significance of differences was determined using ANOVA or two-way ANOVA with Bonferroni correction for post-tests.
Online supplemental material
The supplemental material contains additional information on recruitment of RyR openings in the spark (Eqs. S1–S6); on calcium sensitivity of spontaneous spark frequency (Eq. S7); calcium dependence of calcium sequestration (Eqs. S8 and S9, Table S1, and Fig. S1); on the magnitude of calcium increase at a neighboring release site (Fig. S2); on analysis of calcium release events (Fig. S3); and on the time dependence of calcium release current during the spark (Fig. S4). The online supplemental material is available at http://www.jgp.org/cgi/content/full/jgp.201411281/DC1.
In the cardiac myocyte, formation of calcium waves is reflected in increased calcium spark frequency and in temporal correlation of spark occurrence at neighboring calcium release sites. We conjecture that both aspects of wave formation should be quantitatively related to RyR gating, to organization of RyRs into calcium release sites, and to properties of the cell cytosol that modulate calcium dynamics. To find the relationship between calcium-dependent gating of RyR channels within the calcium release site and properties of calcium release events, we simulated a system of release sites corresponding to a typical line scan, as described in Materials and methods.
Diastolic calcium release under control conditions
Spontaneous diastolic calcium sparks simulated using the aHTG model of the RyR.
Spontaneous calcium sparks were simulated by incorporating our models of RyR gating and calcium spark generation (Zahradníková et al., 2010; Zahradníková and Zahradník, 2012) into a modified stochastic FDF model based on the model of Coombes et al. (2004) with 51 release sites (CRUs) spaced by 2 µm (see Materials and methods). The cytosolic concentration of free Ca2+ was varied in the range of 0.05–0.4 µM. The parameters characterizing the RyR gating model (Table 1) and the CRUs (Table 2) were identical to those shown previously (Zahradníková and Zahradník, 2012) to provide theoretical estimates of RyR open probability and of the frequency of spontaneous calcium sparks that have been shown previously (Zahradníková et al., 2010; Zahradníková and Zahradník, 2012) to be in excellent agreement with the experimental data (Lukyanenko and Gyorke, 1999; Zahradníková et al., 2003). Parameters characterizing Ca2+ diffusion and sequestration (Table 3) were chosen to imitate conditions, at which spontaneous calcium sparks do not propagate into waves, such as in the presence of 500 µM EGTA (DCa = 30 µm2 s−1; Kushmerick and Podolsky, 1969; Allbritton et al., 1992; Baylor and Hollingworth, 1998; Cordeiro et al., 2001). Calcium release flux through a single RyR channel was adjusted to conform to a single-channel current of 0.4 pA (see Materials and methods; Mejía-Alvarez et al., 1999; Janíček et al., 2012). The calcium release flux during the spark was assumed to have constant amplitude of an integer multiple of the unitary flux, corresponding to nO open RyRs (Wang et al., 2004; Janíček et al., 2012). Calcium release flux in all sparks was assumed to last 15 ms (Zahradníková et al., 2007; Janíček et al., 2012).
Visualization of a typical simulation under control conditions and at different cytosolic [Ca2+] is shown in Fig. 1 A. Nonlinearity of calcium uptake at elevated cytosolic calcium levels was taken into account by using a calcium-dependent value of τd (see Fig. S1). The dependence of calcium spark frequency per CRU on free cytosolic [Ca2+] is shown in Fig. 1 B together with the theoretical prediction of spontaneous calcium spark frequency using the RyR gating model (Zahradníková et al., 2010) and with the experimental data taken from Lukyanenko and Gyorke (1999). In the cytosolic [Ca2+] range of 0.05–0.25 µM, the frequency of calcium spark occurrence during the simulations closely followed the experimentally observed and theoretically predicted calcium spark frequency. At 0.4 µM of cytosolic [Ca2+], an increase of the simulated spark frequency above the value theoretically predicted for spontaneous sparks could be observed together with the appearance of calcium waves (Fig. 1 A). This behavior is consistent with experimental observation of aborted calcium waves at 0.4 µM of cytosolic [Ca2+] (Lukyanenko and Gyorke, 1999). However, the presented simulation does not distinguish whether the observed increase of calcium spark frequency above the value predicted for spontaneous sparks is caused by the calcium dependence of RyR gating or by the calcium dependence on the rate of calcium uptake into the SR.
Spontaneous diastolic calcium waves simulated using the aHTG model of the RyR.
In experiments on skinned cardiac myocytes, decreasing the concentration of EGTA from 500 to 100 µM led to initiation and propagation of calcium waves (Lukyanenko and Gyorke, 1999). The conditions of this experiment were mimicked by increasing the calcium diffusion coefficient (Allbritton et al., 1992; Wagner and Keizer, 1994) while keeping the peak calcium concentration during the spark constant, as experimentally observed (Lukyanenko and Gyorke, 1999). Decreased buffering capacity should lead to increased calcium concentration away from the release site, as experimentally observed and theoretically predicted by calculations (Sham, 1997). To achieve this increase, it was necessary to keep the peak calcium concentration constant, which was accomplished by a fold increase of the parameter jCa for an n-fold increase of the parameter DCa. Other parameters of the simulation remained unchanged. Visualization of a typical simulation at different values of n = 1, 2, and 3 is shown in Fig. 2 A. It is apparent that simulated EGTA removal profoundly increases the propensity of calcium sparks to propagate into calcium waves. Formation of calcium waves at low EGTA (n = 3) was significantly higher at low to intermediate cytosolic [Ca2+], as apparent from the calcium dependence of the frequency of triggered calcium sparks (Δϕspark; Fig. 2 B). The average number of open RyRs in spontaneous sparks (i.e., in solitary sparks and sparks that initiated a calcium wave) was not substantially affected by simulated EGTA removal (Fig. 2 C, open symbols). However, wave propagation induced a large and significant increase of the number of open RyRs per CRU, so that the observed average number of open RyRs in triggered sparks increased by up to 80% upon decrease of EGTA (ANOVA; P < 0.001; Fig. 2 C, closed symbols), thus mimicking the increased amplitude of calcium elevation during the calcium wave with respect to the elevation during a spontaneous calcium spark (Lukyanenko and Gyorke, 1999). In further analyses, we have denoted this increase in the average number of open RyRs in triggered sparks, relative to the number predicted for spontaneous sparks under the same conditions, as
It can be concluded that the simulations of diastolic calcium release under control conditions are in a very good agreement with experimental data.
Regulation of calcium wave formation by RyR gating
To separate the effect of RyR gating on the formation of calcium waves from the effect of varying calcium uptake rate, calcium diffusion, or unitary calcium release flux, we have performed a set of simulations in which individual parameters of the model were varied one at a time. The remaining parameters were kept at their standard values as listed in Tables 1–3.
Several reports (see below) have shown that under pathological conditions, calcium waves may occur without a substantial change in the amplitude and duration of spontaneous solitary sparks, and without a change in the experimental conditions. That is, they may occur under conditions when the only changed entity is the gating of the RyR channel. Most usually, the increased occurrence of calcium waves goes hand in hand with increased calcium spark frequency (Terentyev et al., 2006, 2008, 2009; Fernández-Velasco et al., 2009). However, in one particular report, calcium spark frequency in cardiac myocytes from failing hearts was increased but the propensity to form calcium waves was decreased (Song et al., 2005). A mechanistic model of calcium spark formation should be able to reproduce these apparently contradictory results. To find out the relationships between the multi-parameter space of the model and the propensity for initiation and propagation of calcium waves, we have analyzed the relevant model properties.
Formation of a propagating calcium wave requires that the elevation of calcium concentration caused by a calcium spark results in a sufficiently high probability of production of a triggered calcium spark at a quiescent neighboring CRU. For instance, to achieve a 20% probability of wave propagation across five neighboring CRUs, the average probability of spark triggering at these CRUs has to be At the same time, the probability of spontaneous occurrence of a calcium spark at these CRUs has to be adequately low to reproduce the experimentally observed calcium spark frequencies (e.g., for the experimentally observed range of 10–50 sparks per second per 100-µm scanning line, Eq. 10 gives the probability that a 15-ms spark occurs at a single CRU as PCRU = 0.003–0.015). To achieve high triggering probability at low spontaneous calcium spark frequency, the increase of local cytosolic [Ca2+] caused by a spark activated in the vicinity of a CRU has to induce a large increase of the calcium spark frequency. To quantify the calcium sensitivity of spark frequency for different parameter sets, we have expressed it as ∂Caϕspark (see Materials and methods).
Modulation of RyR gating such as by adrenergic stimulation was shown to alter the average number of open RyRs in calcium sparks (Zhou et al., 2009), and such changes in may be explained by altered RyR gating (Zahradníková et al., 2010). Because the calcium release flux is directly proportional to it can be expected to affect the formation of triggered sparks. To quantify the effect of the average number of open RyRs in calcium sparks for different parameter sets, we expressed for each parameter set using Eq. 8.
Additionally, the calcium sensitivity of can be expected to affect wave propensity as well, as increased sensitivity will lead to progressively increasing the number of open RyRs in triggered calcium sparks and consequently to increased probability of further wave propagation. To quantify the calcium sensitivity of the average number of open RyRs in spontaneous calcium sparks, we have expressed it for each parameter set as (see Materials and methods).
Modulation of calcium release by changes in RyR gating parameters.
Calcium sparks were simulated in the range of cytosolic [Ca2+] spanning the range from 0.05 µM, corresponding to intact cells during the diastole (De Young and Scarpa, 1987; Smogorzewski et al., 1998), to 0.4 µM. Changes of parameters related to the RyR activation site—the dissociation constants of Ca2+ and Mg2+ from the activation site (KCa, KMg), the allosteric factors of Ca2+ and Mg2+ (fCa, fMg), and the probability of Mg2+ dissociation from the activation site during a spark (pdis)—are examined in Fig. 3. The remaining parameters were assigned their control values (Tables 1–3). Fig. 3 A shows the results of typical simulations for cytosolic [Ca2+] of 0.05 (top row) and 0.4 µM (bottom row) in control or when the evaluated parameters were changed in the direction that increased calcium spark frequency, i.e., a decrease to 33% for KCa, fCa, and fMg, an increase to 300% for KMg, and an increase to 200% for pdis. Fig. 3 (B and C) shows the values of the frequency of triggered sparks (Δϕspark) and of the increase in the average number of open RyRs in triggered sparks respectively, for control, with increased and decreased parameters in the whole range of diastolic Ca2+ (closed symbols; the meaning of the lines in B and C is explained below). Under control conditions, no wave formation occurs and the values of Δϕspark and are close to 0 in the whole [Ca2+] range, in contrast to data in Fig. 1. This result indicates that under control conditions, the tendency to form waves at elevated cytosolic Ca2+ of 0.4 µM is caused by the decelerated calcium sequestration into the SR and not by RyR gating.
Increase of calcium spark frequency caused by changes in RyR gating resulted in the appearance of calcium waves upon changes in KCa, KMg, fCa, and pdis, but not when fMg was decreased. Whereas modified KCa, KMg, and fCa had the strongest effect on wave formation at low cytosolic [Ca2+] (Fig. 3 B) and resulted in an increased number of open RyRs in triggered sparks (Fig. 3 C), increased pdis led to increased wave formation in the whole cytosolic [Ca2+] range (Fig. 3 B) but had a much smaller effect on the increase of the number of open RyRs in triggered sparks (Fig. 3 C).
Comparison of the calcium dependence of the properties of triggered sparks (Fig. 3, B and C) and the properties of spontaneous sparks (Fig. 3, D–G) shows that spontaneous spark frequency ϕspark (Fig. 3 D) has practically no effect on spark triggering, whereas both Δϕspark and seem to change in proportion to the calcium sensitivity of spontaneous spark frequency, ∂Caϕspark (Fig. 3 E, cf. panels corresponding to changed KCa, KMg, and fCa), and to the average number of open RyRs in spontaneous sparks (Fig. 3 F; cf. the panel corresponding to changed pdis). The calcium sensitivity of (Fig. 3 G), seems to affect the number of open RyRs in triggered sparks but not their frequency, as apparent from the elevation of at [Ca2+] = 0.4 µM after a decrease of the parameters fCa and KCa, but not after an increase in KMg (Fig. 3, cf. C and G).
Changes of the remaining parameters related to spontaneous spark formation—the closed state dissociation constant, KO00 (an inverse measure of RyR tendency to spontaneous opening), the Mg2+ dissociation constant of the inhibition site, KI, and the refractory period trp—are examined in Fig. 4. Fig. 4 A shows the results of typical simulations for cytosolic [Ca2+] of 0.05 (top row) and 0.4 µM (bottom row) when the evaluated parameters were changed in the direction that increased calcium spark frequency, i.e., a decrease to 33% for KO00 and trp and an increase to 300% for KI. Neither of the changes brought about a strong tendency to wave formation (Fig. 4 A). The decrease of KO00 and trp and the increase of KI led to only a moderate increase of the frequency of triggered sparks (Fig. 4 B), despite substantially elevating spontaneous spark frequency (Fig. 4 D). The number of open RyRs in triggered sparks (Fig. 4 B) was increased only slightly as well. The relationships between the properties of triggered (Fig. 4, B and C) and spontaneous sparks (Fig. 4, D–G) paralleled those observed in Fig. 3; i.e., the increase in both the frequency and in the number of open RyRs in triggered sparks was proportional to calcium sensitivity of spontaneous spark frequency (Fig. 4 E) and to the number of open RyRs in spontaneous sparks (Fig. 4 F).
Modulation of calcium release by changes in calcium-handling parameters
As demonstrated by the differences between Figs. 1, 2, and 3, wave propagation is affected not only by changes in RyR gating but also by the calcium-handling properties of the cell, i.e., by the rate of calcium uptake (inversely proportional to the SR uptake time constant τd), by the amplitude of calcium flux through a single open RyR (jCa), and by the apparent calcium diffusion coefficient DCa.
The effect of a threefold change of the parameters DCa, jCa, and τd on calcium sparks and waves is shown in Fig. 5. Panels in this figure are organized as those of Figs. 3 and 4. Increase of all the studied parameters led to the appearance of well-formed waves (Fig. 5 A, top row), which persisted even at increased cytosolic [Ca2+] (Fig. 5 A, bottom row). In contrast to the data presented in Figs. 3 and 4, the highest values of Δϕspark for increased parameters were observed especially at intermediate [Ca2+] (Fig. 5 B). The value of was significantly different from 0 for increased parameters DCa and jCa, but above all for increased τd.
Quantitative relationships characterizing the propensity to form calcium waves
Qualitatively, it is apparent that the highest incidence and persistence of waves occurred under conditions that increased frequency of triggered sparks, Δϕspark. Increase of spark frequency per se, either by a change of RyR gating (parameter fMg in Fig. 3 A and parameters KO00, KI, and trp in Fig. 4 A) or by increased cytosolic [Ca2+] (bottom panels in Figs. 3 A and 4 A) was not sufficient for the appearance of calcium waves.
To quantify the relationships between the variables describing the properties of spontaneous calcium sparks (ϕspark, ∂Caϕspark, , input variables characterizing calcium dynamics of the cell (DCa, jCa, and τd), and the variables characterizing the triggered calcium release events (Δϕspark, , we have performed multiple linear regression analysis.
The frequency of triggered sparks, Δϕspark, was dependent on the product ∂Caϕspark and on the weighted average of the squares of variables characterizing calcium dynamics, DCa, jCa, and τd, resulting in the equation(13)which described the data with coefficient of determination (R2) = 0.92 and coefficient of variation (ĉv) = 0.53 (solid lines in Figs. 3 B, 4 B, and 5 B). The proportionality of Δϕspark to ∂Caϕspark is caused by the commensurability between and the local [Ca2+] elevation at the site of the triggered spark. The lack of a direct dependence of Δϕspark on spontaneous spark frequency can be explained by the implicit proportionality between ϕspark and ∂Caϕspark. The quadratic dependence of Δϕspark on the calcium-handling parameters (diffusion coefficient DCa, unitary calcium release flux jCa, and time constant of calcium sequestration τd) indicates nonlinear character of the spark-triggering process. The parameters A1–A3 describing the relative effect of the variables DCa, jCa, and τd are summarized in Table 4. The results are shown as solid lines in Figs. 3 B, 4 B, and 5 B.
The increase of the average number of open RyRs in triggered sparks, was most strongly dependent on the sensitivity of to changes in cytosolic calcium It was also dependent on the product ∂Caϕspark and on a weighted function of the variables characterizing calcium dynamics, DCa, jCa, and τd. The best description of the data was achieved by assuming a quadratic dependence of on the diffusion coefficient DCa, and fourth-order dependence on the unitary calcium release flux jCa and on the time constant of calcium sequestration τd, resulting in the equation(14)which described the data with R2 = 0.88 and ĉv = 0.46 (solid lines in Figs. 3, 4, and 5 D). The higher order dependence of on the calcium-handling parameters again indicates nonlinear character of the spark-triggering process. The parameter A0 describing the effect of the variable and the parameters A1–A3 describing the relative effect of the variables DCa, jCa, and τd are summarized in Table 5. The results are shown as solid lines in Figs. 3 C, 4 C, and 5 C.
Despite the simplicity of the equations and the large variety of the observed relationships between RyR gating parameters, cytosolic [Ca2+], and the properties of triggered calcium sparks, the description of the calcium dependence of Δϕspark and by Eqs. 13 and 14 is very good.
The relative effect of individual properties of spontaneous sparks (∂Caϕspark, and ) and of individual calcium-handling parameters (DCa, jCa, and τd) on triggered calcium sparks was analyzed using the predictions of Eqs. 13 and 14 for their individual input variables, whereas the effect of the remaining variables was suppressed. To suppress the effect of ∂Caϕspark and their values were kept equal to those in the absence of cytosolic Ca2+. To suppress the effect of DCa, jCa, and τd, the parameters of Eqs. 13 and 14—A0, A1, A2, and A3, respectively—were set to zero. The impact of individual properties of spontaneous sparks on the properties of triggered sparks is shown in Fig. 6. Fig. 6 A indicates that the frequency of triggered calcium sparks (Δϕspark) is controlled to a large extent by the calcium sensitivity of calcium spark frequency (∂Caϕspark). The number of open RyRs in spontaneous sparks contributes to the frequency of triggered sparks most strongly upon changes in the parameter pdis, when is the only changed property of spontaneous sparks, and significantly also upon changes in the parameters fCa and KMg in the low cytosolic [Ca2+] range. The relative effect of ∂Caϕspark and on the increase of the number of open RyRs in triggered sparks Fig. 6 B) is less prominent, as the latter is dominated by the effect of especially at high cytosolic [Ca2+]. At low cytosolic calcium, the effect of ∂Caϕspark may contribute up to 50% to under conditions when this parameter is responsible for wave formation (i.e., upon changes in parameters controlling calcium binding to the activation site). A high contribution of to was observed only upon increase in the probability of Mg2+ dissociation during the spark, which affects but does not change the remaining properties of spontaneous sparks.
Fig. 7 depicts the impact of individual calcium-handling parameters on the properties of triggered sparks. Fig. 7 A shows that the frequency of triggered sparks is affected the most strongly by the time constant of calcium reuptake into the SR, τd. This effect is predominant under most conditions, but when DCa or jCa is increased, its effect becomes similar in size or even larger than that of τd. The relative contribution of the calcium-handling parameters DCa, jCa, and τd is not significantly dependent on cytosolic calcium concentration. The impact of individual calcium-handling parameters on the increase in the number of open RyRs in triggered sparks is shown in Fig. 7 B. Upon changes in parameters of the RyR gating model (fCa, pdis), the contribution of DCa to the increase of the number of open RyRs in triggered sparks Fig. 7 B) was larger than that of jCa and τd. Increase in any of the calcium-handling parameters strongly increased their effect on as expected from their higher order dependence (Eq. 14).
To examine whether the formation of waves is solely caused by the appearance of triggered sparks, or whether sparks may become organized into waves without an increase in their total frequency, we have tested the presence of correlation between the frequency of triggered sparks (Δϕspark) and the average number of sparks in wave events Fig. 8 A). If organization of sparks into waves without increase in spark frequency is possible, we would expect that at least under some conditions an increase in without a corresponding increase of Δϕspark should occur. For clarity, different types of parameter changes were grouped together and color-coded: LA (short for low activity; black) represents control conditions and parameter changes that decreased the frequency of spontaneous calcium sparks or the spatial spread of Ca2+; RyR-Act (red) represents parameter changes that increased spontaneous calcium spark frequency by modulating the RyR activation site; RyR (green) represents parameter changes that increased spontaneous calcium spark frequency by modulating other properties of RyR gating; and CaH (short for Ca handling; blue) represents parameter changes that increased the spatial spread of Ca2+. For groups RyR-Act and CaH, there was a highly significant correlation (Pearson’s correlation coefficient [RP] = 0.90 and 0.89, respectively; P < 0.001) between Δϕspark and the slopes of which were not significantly different (4.56 ± 0.47 and 4.93 ± 0.71 s, respectively). The correlation between Δϕspark and in the RyR group was marginally significant (RP = 0.47; P = 0.08), and its slope was much shallower (1.22 ± 0.16 s; green line in Fig. 8 A). We interpret the decreased slope as indicating that under conditions of the highest spontaneous spark frequency (decreased parameter trp, [Ca2+] = 0.25 and 0.4 µm; labels in Fig. 8 A), a fraction of spontaneous sparks was erroneously categorized as triggered sparks, as a result of their being separated by <50 ms from the spark occurring at a neighboring location. Under conditions of low spontaneous spark frequency (LA), the data points were clustered near the origin of the graph (Δϕspark = 0; ), and no correlation was observed (RP = 0.44; P = 0.18). Overall correlation between Δϕspark and for all data were also highly significant (RP = 0.91; P < 0.001; slope, 4.88 ± 0.28 s). Thus, the parameter Δϕspark is suitable for quantification of the propensity to form calcium waves, and no additional organization of sparks into waves was observed in the simulations.
The average number of sparks in the wave was relatively low for conditions with the highest propensity for calcium waves). This is in part caused by the geometric distribution of nSpark of individual waves (not depicted). However, even under conditions supporting the highest wave formation (decreased fCa or increased τd), waves formed by 20 or more sparks were rarely observed (approximately one to two per 10 s), and only one wave with 30 or more sparks occurred in the whole course of simulations (not depicted). To obtain further insight into the mechanisms of wave termination, the dependence of calcium wave termination on RyR gating parameters, variables characterizing calcium dynamics of the cell, and calcium concentration was examined by two-way ANOVA. Simulations were categorized into four groups as in the previous analysis. The results are shown in Fig. 8 (B–E). Termination of waves by propagation failure and by refractoriness depended on both the type of parameter change and calcium concentration, whereas termination by collision of waves depended only on the type of parameter change (two-way ANOVA; P < 0.05). There was no interaction between the effect of cytosolic Ca2+ concentration and the effect of the remaining simulation parameters on the mechanism of calcium wave termination. The fraction of waves terminated by propagation failure decreased (RP = −0.99; P < 0.002) and the fraction of waves terminated by refractoriness increased (RP = 0.99; P < 0.002) with cytosolic [Ca2+]. Propagation failure was the largest in the LA group; it was significantly smaller in the RyR and CaH groups; and it was the smallest in the RyR-Act group. Termination by refractoriness was significantly larger in the RyR-Act group than in the remaining groups. Termination by wave collision was the largest for the RyR-Act group; in the CaH group, it was significantly larger than in the LA group but not significantly different from the RyR group; and it was not significantly different in the RyR and LA groups. Collectively, the data in Fig. 8 (B–E) indicate that under conditions supporting wave formation, waves are terminated by refractoriness and/or wave collision because of the high frequency of calcium release events.
Suppression of calcium waves by controlling the spontaneous RyR opening transition
We have shown that impairment of RyR gating properties may lead to calcium waves as a result of increased calcium sensitivity of spark frequency (∂Caϕspark) or increased number of open RyRs in sparks This begs the question of whether it is possible to prevent wave formation induced by diverse RyR gating changes (i.e., by increased ∂Caϕspark vs. by increased ) by a common treatment at the level of RyR gating. We have examined the effect of a gating parameter that by itself has a weak effect on formation of calcium waves: the spontaneous RyR opening transition (increase of the RyR gating parameter KO00). The effect of increased KO00, which impeded spontaneous RyR opening, was tested. Calcium waves were induced by increasing the calcium affinity of the activation site (decrease of the RyR gating parameter KCa), decreasing the magnesium affinity of the activation site (increase of the RyR gating parameter KMg), increasing the positive allosteric effect of Ca2+ binding to the RyR activation site (decrease of the RyR gating parameter fCa), or increasing the probability of Mg2+ unbinding (increase of the RyR gating parameter pdis). The effect of impeding spontaneous opening of RyR on calcium waves that have been induced by the above described ways of RyR activation site modulation is illustrated in Fig. 9. Typical simulations under control conditions, after inducing wave formation by decreased fCa, and after combining the effect of decreased fCa with the effect of increased KO00, is shown in Fig. 9 (A–C). The formation of waves arising because of a threefold decrease of fCa (cf. Fig. 9, A and B) can be almost fully reversed by a threefold increase in KO00 (Fig. 9 C). Statistical analysis for all tested RyR gating parameters is shown in Fig. 9 D.
The ability of changes in KO00 to suppress wave formation induced by diverse RyR gating changes suggests that the spontaneous opening transition may be a useful drug target to treat calcium-dependent arrhythmias caused by malfunction of RyR gating arising from different causes.
In this work, we present a framework for stochastic simulation of calcium release events based on a model of RyR gating that has been shown previously to describe RyR gating in bilayers as well as the properties of calcium sparks in cardiac myocytes (Zahradníková et al., 2010; Janíček et al., 2012; Zahradníková and Zahradník, 2012), and at physiological levels of calcium release flux (∼0.4 pA per open RyR; Mejía-Alvarez et al., 1999; Janíček et al., 2012).
Our data show that the propensity to form calcium waves does not have a direct relationship to spontaneous spark frequency: it may increase with increased spontaneous spark frequency (e.g., by changes in gating that enhance calcium binding to the RyR activation site: by a decrease in KCa or fCa, or an increase in KMg at 0.05 µM [Ca2+]; Fig. 3); it may increase without a change of spontaneous spark frequency (e.g., by changes in gating that enhance magnesium unbinding from the RyR activation site by an increase in pdis at 0.05 µM [Ca2+]; Fig. 3); it may remain unchanged after an increase in the spontaneous spark frequency (e.g., a decrease in fMg at 0.25 µM [Ca2+]; Fig. 3), but it may also decrease after an increase in the spontaneous spark frequency (e.g., an increase of [Ca2+] at fCa = 0.015; Fig. 3), which may explain experimental results that linked increased spark frequency to decreased wave propensity in compensated left ventricular hypertrophy (Song et al., 2005) and the decreased sensitivity of RyRs to activators in some arrhythmogenic RyR mutations (Thomas et al., 2004). The presented simulations suggest that changes in RyR gating may have profound and diverse effects on formation and propagation of waves.
Based on theoretical considerations, we have speculated that diastolic calcium wave formation upon RyR2 malfunction may be a consequence of elevated calcium sensitivity of calcium spark frequency and increased average number of open RyRs in spontaneous sparks. We have confirmed this proposal using computer simulation of calcium sparks. The propensity to form calcium waves was further modulated by the calcium-handling properties of the cell, and the number of open RyRs in triggered sparks was modulated by the calcium sensitivity of the average number of open RyRs in spontaneous sparks. The combined effect of the above factors resulted in a high propensity to form calcium waves at diastolic (0.05–0.1 µM) cytosolic [Ca2+] when the parameters affecting the calcium-binding properties of the RyR activation site were changed in the direction that increased RyR open probability and frequency of spontaneous calcium sparks (KCa, KMg, fCa). When calcium-binding properties of the RyR activation site were changed in the direction that increased the number of open RyRs in spontaneous calcium sparks (pdis), a high propensity to form calcium waves could be observed in the whole examined range (0.05–0.4 µM) of cytosolic [Ca2+].
Calcium wave formation at low frequency of spontaneous calcium sparks
To our knowledge, this is the first report that uses the same theoretical framework for quantitative description of the formation of spontaneous calcium sparks and the initiation and propagation of calcium waves. Previous mathematical models often used a predefined threshold of cytosolic [Ca2+] to start a calcium wave (Dupont et al., 1996; Keizer et al., 1998; Lukyanenko et al., 1999; Subramanian et al., 2001; Okada et al., 2005). Because activation of calcium release is a stochastic process (Cheng et al., 1993; Cannell et al., 1994; Poláková et al., 2008; Janíček et al., 2012), efforts have been made to incorporate stochasticity into these models by introducing a simple stochastic channel-gating scheme (Keizer and Smith, 1998), a stochastic waiting-time distribution (Izu et al., 2001, 2006), or a stochastic activation threshold (Coombes and Timofeeva, 2003; Coombes et al., 2004).
However, previous works mostly concentrated on the calcium dynamics during calcium sparks and waves and therefore used phenomenological models of calcium release site activation that have not been tested for their conformity with the properties of single RyR channels. Moreover, these models often had to use separate sets of conditions for simulating spontaneous calcium sparks and propagating calcium waves, as initiation of waves required simultaneous activation of several neighboring calcium release sites to ensure successful propagation of the wave (Izu et al., 2001, 2006), whereas under diastolic conditions such a simultaneous activation is highly unlikely. For instance, at a calcium spark frequency of 25 s−1 (100 µm)−1 and for a spark duration of 15 ms, the probability that two CRUs out of the nine closest neighbors of an active CRU in the transversal plane of a cardiac myocyte will simultaneously activate is P = 5.6 × 10−5. Under many pathological conditions, the frequency of diastolic calcium sparks is only moderately elevated (Cheng et al., 1993; Maier et al., 2003; Kubalova et al., 2005; Guo et al., 2006; Terentyev et al., 2006, 2007; Fernández-Velasco et al., 2009), but the occurrence of spontaneous calcium release was shown (Cheng et al., 1993; Terentyev et al., 2006, 2007; Fernández-Velasco et al., 2009) or supposed (Yano et al., 2006) to underlie the studied pathologies (Yano et al., 2006). Paradoxically, even a decrease in the propensity to form calcium waves has been observed in parallel to elevated calcium spark frequency (Song et al., 2005). These observations suggest that the presence of several adjacent active release sites is not necessary for calcium wave initiation.
Here, we made use of the findings that if calcium sparks are allowed to occur only at the points of a discrete, regular temporal lattice and are assumed to have the form of a rectangular pulse, the spatiotemporal evolution of each calcium spark can be calculated analytically, thus eliminating the need to solve the system of reaction–diffusion equations (Izu et al., 2001; Coombes et al., 2004). This principle enabled the analysis of the relationship between RyR activity, generation of spontaneous sparks, and generation of calcium waves using a mechanistic model of RyR gating and calcium spark formation (Zahradníková et al., 2010; Zahradníková and Zahradník, 2012) and realistic estimates of the composition of CRU (Zahradníková and Zahradník, 2012). The linear formulation of the model, although a simplification, enabled us to separate the effect of RyR gating and those of cytosolic calcium handling and examine them in isolation.
We have found that indeed, calcium waves may be initiated at a low calcium spark frequency if the sensitivity of the changes in calcium spark frequency to Ca2+ concentration is sufficiently elevated. The appearance of waves was consistent with experimental data (Lukyanenko and Gyorke, 1999) at reasonable values of single-channel calcium release flux and apparent calcium diffusion coefficient. We attribute the ability of our model to generate calcium waves at low calcium spark frequency and with low calcium release fluxes to two factors: first, the use of an exact relationship for the open probability of RyRs that is able to convey the subtle differences of the cytosolic calcium dependence of RyR open probability from the Hill equation; and, second, quantitative incorporation of the presence of many RyRs at the calcium release sites. This approach permitted a realistic description of calcium wave initiation that did not require unrealistically high values of the refractory period, as reported previously for a Hill equation–based model (Izu et al., 2001).
Relationships between spontaneous calcium sparks and calcium waves
The effects of spontaneous calcium spark properties.
Changes of different parameters of RyR gating had disparate effects on the calcium dependence of calcium wave formation (cf. Figs. 3 and 4). Analysis of the relationships between the theoretical properties of spontaneous calcium sparks and variables characterizing calcium dynamics of the cell on one hand and the production of triggered calcium sparks on the other hand revealed that the whole range of observed effects is brought about by a common mechanism. The production of triggered calcium sparks, characterized by the variable Δϕspark, as well as the production of additional RyR openings in triggered calcium sparks, characterized by the variable could be unequivocally described by equations involving solely the properties of spontaneous calcium sparks (∂Caϕspark, ) and parameters characterizing calcium dynamics of the cell (DCa, jCa, and τd). That is, quantifying the production of triggered sparks required only the knowledge of empirically measurable properties of spontaneous calcium sparks without necessitating the knowledge of the RyR gating mechanism. The analysis has shown the crucial role of both the calcium sensitivity of spontaneous spark frequency (∂Caϕspark) and of the number of open RyRs in spontaneous sparks for calcium wave formation.
RyR gating and wave termination.
Statistical analysis revealed that the frequency of triggered calcium sparks (Δϕspark) unequivocally determined the average length of waves . That is, we have not observed any tendency of sparks to organize themselves into waves without increase in the total number of sparks.
Wave termination was dependent on both cytosolic Ca2+ concentration and on the mechanism of triggered spark production. The probability of termination by propagation failure decreased and the probability of termination by refractoriness increased with cytosolic Ca2+ in a complementary fashion, whereas the probability of wave termination by collision of two waves changed to a much lesser extent. The probability of wave termination by propagation failure was decreased and probability of wave termination by refractoriness was increased by all conditions that augmented wave formation, but most of all when RyR open probability was increased by modulation of the RyR activation site.
Gating model limitations
The gating model (Zahradníková et al., 2010) does not have any parameters directly related to control of RyR activity by luminal calcium. It was not possible to construct a quantitative model of luminal calcium regulation because of the scarcity and partial inconsistency of the experimental data (Györke and Györke, 1998; Györke and Terentyev, 2008; Qin et al., 2008; Tencerová et al., 2012; Chen et al., 2013; Gaburjakova and Gaburjakova, 2014). As we suggested previously (Zahradníková et al., 2010; Tencerová et al., 2012), luminal calcium may be expected to influence RyR gating by affecting one or more of the parameters fCa, fMg, and KO00. Out of these parameters, only a change in the allosteric factor fCa increases the propensity to form calcium waves.
No matter whether the refractoriness of calcium release sites is caused by inactivation of RyRs (Sham et al., 1998) or by luminal Ca-dependent deactivation (Terentyev et al., 2003; Zima et al., 2008), it is ultimately caused by changes in RyR gating; therefore, it should be a stochastic variable. We have recently shown (Janíček et al., 2012) that, albeit with a very low probability, a calcium release site may not become refractory after a calcium release event. However, in the absence of theoretical understanding of the principles governing RyR refractoriness, it is at present not possible to determine the probability distribution of refractoriness in individual release sites. Therefore, in the present model, the refractory period was modeled as a constant. This was motivated by the steep, almost step-like time dependence of the probability of calcium spark activation during recovery from refractoriness, observed experimentally (Sobie et al., 2005). Incorporation of luminal regulation of the refractory period will be the subject of our further study.
Steady-state approximation of Mg2+ dissociation from RyR activation sites.
In the presented model, the activation sites of closed RyR channels are assumed to be in rapid equilibrium with the cytosolic Ca2+, which dynamically changes in response to calcium spark formation. This is an oversimplification because of the large extent of saturation of the activation sites by Mg2+ under diastolic conditions and the finite rate of Mg2+ dissociation from the activation sites (Zahradníková et al., 2010; Janíček et al., 2012; Zahradníková and Zahradník, 2012). The bias that might have been introduced into the simulations by this assumption was minimized by keeping the simulation time step at a sufficiently large value (Δt = 3 ms).
Calcium signaling limitations
Time course of RyR single-channel calcium current.
The amplitude of single-channel current through an open RyR, iCa, was assumed constant throughout the spark. However, it can be assumed that iCa will decrease in parallel with the decrease of the calcium content in the SR cisterna. We have tested the effect of time-dependent calcium release flux on calcium spark activation. The procedure is described in the supplemental text. Fig. S4 depicts two sets of simulations. In the first simulation, iCa stayed constant throughout the spark, whereas in the second simulation, it decreased linearly. The time integral of iCa was 6 fC in both simulations. It is apparent that both simulations produce a virtually identical pattern of calcium sparks under control conditions as well as in the presence of wave formation. It can be concluded that the simplification of constant iCa does not affect the conclusions of this work.
The use of the simplified one-dimensional spatial arrangement greatly sped up the calculations and thus allowed us to use complex equations to calculate the probabilities necessary for determination of the occurrence of calcium sparks and of their calcium release flux. We consider this model relevant, because calcium waves are most often experimentally observed in the longitudinal direction, and because the calcium flux in the model was carefully adjusted to conform to the three-dimensional situation, as described previously for the adjustment of two-dimensional models (Izu et al., 2001). Extension of the model to two- and three-dimensional situations will be the topic of our further work.
The presented model does not involve binding of Ca2+ ions to mobile and immobile buffers, which are known to influence cell calcium dynamics. Extension of the model to include such buffers explicitly would require a numerical approach to solve partial differential equations describing calcium diffusion, different from the analytical approach used here. In our simulations, we have assumed that cytosolic calcium buffering can be expressed as a decrease in the apparent diffusion coefficient DCa (Keizer et al., 1998; Coombes et al., 2004). This approximation leads to additivity of calcium concentration increases from different sources, a condition that is not fulfilled if the cytoplasmic calcium buffers become saturated (Maeda et al., 1999), such as when several neighboring release sites are active simultaneously (Izu et al., 2001, 2006). In our hands, a change of buffering capacity and a corresponding change of the apparent diffusion coefficient had opposite effects on the calcium concentration at the neighboring calcium release site, which could be compensated by keeping the peak calcium concentration constant (see section, Spontaneous diastolic calcium waves simulated using the aHTG model of the RyR).
Additionally, the model of calcium diffusion (Eq. 9) assumes that the rate of calcium uptake to the SR is directly proportional to cytosolic calcium concentration. This assumption holds only at low cytosolic [Ca2+], as the cytosolic calcium-binding site of SERCA has submicromolar KCa. Taking into account the time course of calcium increase at the neighboring CRU in response to spark activation (Fig. S2) and the calcium dependence of τd (Fig. S1), it is apparent that under control conditions, the SERCA pump operates in the linear range during wave activation. However, τd is significantly increased at elevated [Ca2+] (Fig. S1). We expect that during initiation of a diastolic calcium wave, the linear approximation is valid because of the relatively low calcium spark frequency. However, as the calcium wave starts to propagate in three dimensions, the released calcium may saturate the cytosolic calcium buffers and cytosolic SERCA calcium-binding sites, and therefore strong deviations from linearity may occur.
Triggered activity caused by spontaneous calcium release is the mechanism suspected to underlie the pro-arrhythmic effects of RyR dysregulation or RyR mutations underlying several genetic or acquired cardiac diseases (Durham et al., 2007). At the level of cardiac myocytes, spontaneous release is observed in the form of calcium waves (Cheng et al., 1996; Lukyanenko et al., 1999; Fernández-Velasco et al., 2009). The present study provides a simple platform to explain the effect of physiological and pathophysiological modulations affecting calcium release on the propensity to propagate calcium waves. Thus, it offers a simple means of how to test whether a proposed mechanism is consistent with the experimental observations.
The study makes important predictions regarding possible mechanisms of alterations of RyR gating and SR calcium homeostasis that may lead to diastolic calcium wave formation, i.e., spontaneous calcium release. Most importantly, our simulations show that there is no direct correspondence between the effect of RyR modulation/dysregulation on the frequency of spontaneous calcium sparks and on the propensity for calcium waves. Changes that manifest as an increase in the apparent calcium sensitivity of RyRs but do so without affecting the RyR activation site properties do not change the propensity to form calcium waves or increase it only marginally. These results are in line with the observations that increasing the RyR open probability is not sufficient to induce arrhythmias (Venetucci et al., 2007). Our simulations have disclosed that the rate of SR uptake (proportional to 1/τd) has a very strong effect on calcium wave propagation and thus are in excellent agreement with the recently observed strong inhibitory effect of phospholamban knockout on wave propagation (Bai et al., 2013).
Cytosolic versus luminal regulation of the propensity to form calcium waves.
There is a relative lack of information on the molecular mechanisms by which arrhythmogenic RyR mutations lead to calcium release. At the cellular level, many mutations have been shown to decrease the threshold luminal calcium concentration for spontaneous store overload–induced calcium release in transfected cells (Jiang et al., 2004). However, as exemplified by the G230C mutation at the level of single channels, a change in both cytosolic and luminal calcium sensitivity (Liu et al., 2013) or only in cytosolic calcium sensitivity (Meli et al., 2011) has been observed. It has been questioned whether a subtle change in calcium sensitivity can produce sufficient spontaneous calcium release to be arrhythmogenic (Liu et al., 2013). Our results show that calcium waves can be induced by RyR gating changes that increase the apparent cytosolic calcium sensitivity of the channel to a similar extent as observed by Meli et al. (2011) by one of three mechanisms: changing calcium or magnesium affinity of the RyR activation site, or allosterically. The most prominent effect occurs when the allosteric coupling between Ca2+ binding to the RyR activation site and channel opening is increased, which can be described by a decrease of the parameter fCa in the aHTG model of RyR gating. Several arrhythmogenic RyR mutations (P2328S, Q4201R, V4653F) were shown to exhibit two to three times decreased Mg2+ sensitivity (i.e., increased KMg) of the RyR activation site (Lehnart et al., 2004). Our simulations show that such a change would be sufficient to induce formation of calcium waves.
Although luminal calcium has a dramatic effect on RyR open probability (Laver, 2007b; Laver and Honen, 2008; Tencerová et al., 2012) and arrhythmogenic mutations may affect luminal calcium sensitivity (Liu et al., 2013), we propose that for increased wave propensity, the effect of luminal calcium has to transpire through modulation of cytosolic calcium sensitivity, as calcium waves propagate by diffusion of cytosolic calcium (Lukyanenko et al., 1999; Swietach et al., 2010). The presented model can explain the common scenario of store overload–induced calcium release (Jiang et al., 2004) by a luminal [Ca2+]-induced decrease of fCa proposed previously (Zahradníková et al., 2010; Tencerová et al., 2012). Alternatively, the change of luminal calcium sensitivity of RyRs observed in cardiovascular diseases may be caused by an allosteric effect of changes in the properties of the RyR activation site on luminal calcium sensitivity (Tencerová et al., 2012). It is questionable whether it is possible to distinguish whether the primary cause of RyR dysregulation is modulation of luminal or of cytosolic calcium binding, as allosteric effects of ion binding tend to modulate the whole set of state transitions of an ion channel.
Moreover, in the physiological range, calcium release flux is supposed to be proportional to the calcium gradient between SR and cytosol (Stern et al., 1999; Shannon et al., 2004). Literature data suggest a steep dependence of calcium leak on calcium load at elevated luminal Ca and a flat dependence at decreased luminal Ca (Shannon et al., 2002). Song et al. (1997) showed that spark amplitude is directly proportional to SR calcium load, and that spark frequency is independent of SR calcium load in the range of 30–100% of control. This is in line with the results of Chen et al. (2013) and Tencerová et al. (2012) who showed that under physiological conditions (in the presence of Mg2+), RyR open probability is independent of luminal calcium concentration in the range of 0–1 mM and at diastolic cytosolic calcium levels (0.05–0.1 µM). Thus, simulation of the effect of changing calcium load by change in the single-channel current iCa may be adequate. Our results show that decreased single-channel calcium release flux caused by decreased SR calcium load does not affect spark frequency, only spark amplitude, whereas increased SR calcium load, by increasing the single-channel calcium release flux, increases both spark frequency and amplitude, and leads to calcium waves even in the absence of modulation of RyR open probability by luminal Ca2+ (Fig. 5), in agreement with the results of Venetucci et al. (2007).
Reversal of susceptibility to calcium waves by inhibition of spontaneous RyR opening.
Direct modulation of calcium and magnesium interaction with the RyR activation site (changes in the binding properties of the activation site, defined by KCa, KMg) and modulation of RyR activation by Ca2+ in indirect fashion (changes in the allosteric coupling between Ca2+ binding and RyR opening, defined by fCa) have a very similar effect on the propensity to form calcium waves; however, the functional elements responsible for direct and indirect modulation may be supposed to reside on different parts of the protein. For instance, changes in calcium affinity of the activation site detected after mutation of the C-terminal amino acid E3987 (Li and Chen, 2001) are consistent with a change in the calcium-binding properties of the activation site (Zahradník et al., 2005); a similar change in RyR calcium sensitivity observed after mutation of the N-terminal amino acid G230 (Meli et al., 2011; Liu et al., 2013) was proposed to change the interaction between the A and B domains of the N-terminal RyR region (Liu et al., 2013); and mutation of the central amino acid P2328 with a similar effect on RyR gating (Lehnart et al., 2004) was proposed to alter the interaction of RyR with ATP (Blayney et al., 2013), which is compatible with allosteric effects on the occupancy of the RyR activation site (Tencerová et al., 2012). Additionally, a different type of modulation of the same site may produce a different pattern of spark activation and wave formation, such as a change in the Mg2+ dissociation constant (KMg) versus its unbinding rate (pdis) from the RyR activation site. It is therefore of interest to find out whether such mechanistically or functionally disparate changes of RyR activity may be remedied by a common approach. We have tested the effect of inhibiting the spontaneous opening transition of the RyR on calcium waves evoked by four RyR parameter changes differing both in the site of action and in the functional consequences (KCa, KMg, fCa, and pdis). The results are reassuring, as inhibition of spontaneous RyR opening deeply and significantly suppressed the formation of calcium waves in all studied cases.
In brief, this study shows that changes in RyR gating may lead to the appearance of calcium waves at diastolic cytosolic [Ca2+] at physiological levels of single-channel RyR calcium current and with size and spacing of calcium release sites corresponding to experimental findings. This result was obtained by incorporation of an RyR gating model that is in quantitative agreement with data on single-channel RyR open probability, spontaneous calcium spark frequency, and the extent of activation of calcium release sites during cardiac excitation coupling. Analysis of initiation and propagation of waves has shown that for generation of waves at diastolic calcium concentrations, either the sensitivity of calcium spark frequency to changes in cytosolic [Ca2+] or the average number of open RyRs in spontaneous calcium sparks has to be increased. Gating changes fulfilling these conditions occur upon changes in properties of the RyR activation site (KCa, KMg, fCa, pdis) and upon changes in cytosolic [Mg2+]. Besides RyR gating, calcium waves were shown to depend heavily on calcium release flux through single RyRs and on the rate of SR calcium uptake, and less strongly also on the apparent calcium diffusion coefficient. The whole range of observed effects was found to be brought about by a common mechanism that depended on the interaction between ∂Caϕspark, and parameters characterizing calcium dynamics of the cell (DCa, jCa, and τd).
The input of Ivan Zahradník at the early stages of the work and his help in deriving Eqs. 4–6 is gratefully acknowledged. We also thank Steve Coombes and Yulia Timofeeva for providing us with the source code of their stochastic FDF model.
This work was supported by the European Union (contracts LSHM-CT-2005–018802/CONTICA and LSHM-CT-2005-018833/EUGeneHeart), by grants APVV-0721-10 and VEGA 2/0148/14, and by the European Science Foundation program FuncDyn. The funding sources had no involvement in study design; in the collection, analysis, and interpretation of data; in the writing of the report; and in the decision to submit the article for publication.
The authors declare no competing financial interests.
Author contributions: P. Petrovič, I. Valent, E. Cocherová, and A. Zahradníková wrote the MATLAB simulation program; I. Valent and A. Zahradníková derived the theoretical equations; P. Petrovič, I. Valent, E. Cocherová, J. Pavelková, and A. Zahradníková performed simulations; E. Cocherová and A. Zahradníková wrote the MATLAB analysis programs; P. Petrovič, E. Cocherová, J. Pavelková, and A. Zahradníková analyzed the simulations; and I. Valent and A. Zahradníková wrote the paper. All authors edited the paper and acknowledged the final form of the paper.
Eduardo Rios served as editor.
- Abbreviations used in this paper:
- aHTG model
- allosteric homotetrameric gating model
- calcium release unit
- sarco/endoplasmic reticulum calcium ATPase
- Submitted: 27 August 2014
- Accepted: 10 April 2015
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/).