## Abstract

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 Ca^{2+} 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 Ca^{2+} 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.

## INTRODUCTION

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 Mg^{2+} 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 “Ca^{2+}-induced Ca^{2+} 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

### Theory

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 *P _{O}* can be expressed as:

*K*and

_{Ca}*K*are the Ca

_{Mg}^{2+}and Mg

^{2+}dissociation constant of the ion-free closed-state C00;

*f*and

_{Ca}*f*are the allosteric factors coupling Ca

_{Mg}^{2+}and Mg

^{2+}binding to channel opening;

*K*

_{O}_{00}is the equilibrium constant of the O00–C00 transition (inversely proportional to the spontaneous opening tendency of the RyR);

*K*,

_{OL}*K*, and

_{CL}*K*are the dissociation constants of the L-mode states OL and CL and the I-mode state I;

_{CI}*K*is the Mg

_{I}^{2+}dissociation constant of the divalent inhibitory site; and [Ca

^{2+}] and [Mg

^{2+}] 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):*ϕ _{spark}* is the frequency of spontaneous sparks per CRU,

*t*denotes the duration of release,

_{rd}*t*is the period during which the CRU remains refractory after calcium release,

_{rp}*τ*represents the mean RyR open time, and

_{o}*n*is the median number of RyR channels in the CRU (Zahradníková and Zahradník, 2012).

_{RyR}The average period between spontaneous sparks at a single CRU can be expressed as the sum *t _{rd}* +

*t*+

_{rp}*t*, where the idle time

_{it}*t*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

_{it}*τ*of a single RyR and with the median number of RyRs in the CRU,

_{C}*n*(Kunze et al., 1985; Zahradníková et al., 2010; Zahradníková and Zahradník, 2012):

_{RyR}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 *n _{O}* RyR channels are open during a spark,

*Pn*, can be determined using the binomial equation:

_{O}*p*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

_{O}*p*can be calculated from the open probabilities of RyRs with a defined Mg occupancy and from the probability,

_{O}*p*, that a Mg

_{dis}^{2+}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

*p*calculation is given in the supplemental text (Eq. S1). Besides its use in Eq. 7, the probability

_{O}*p*allows for determination of the average number of RyRs that open during a spontaneous spark,

_{O}_{Ca}

*ϕ*) and of the average number of open RyRs in spontaneous calcium sparks

_{Spark}^{2+}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 Ca^{2+} concentration in the vicinity of nearby CRUs. As the open probability of RyRs strongly depends on [Ca^{2+}], 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 Ca^{2+} ions of duration *t _{rd}* and amplitude

*J*. After termination of release, the CRU becomes refractory for the period

_{Ca}*t*. Both

_{rp}*t*and

_{rd}*t*are restricted to integer multiples of Δ

_{rp}*t*. The amplitude

*J*represents the physical flux of Ca

_{Ca}^{2+}(in mol m

^{−2}s

^{−1}). Calcium is released uniformly to both directions. The flux

*J*is given by the unitary flux

_{Ca}*j*, 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):

_{Ca}In Eq. 9, *n _{CRU}* is the number of CRUs;

*m*is the number of calcium sparks occurring at the

_{i}*i*th CRU;

*j*th calcium spark of the

*i*th CRU occurred;

*u*denotes Ca

^{2+}concentration at time

*t*and location

**r**;

*j*th calcium spark at the

*i*th CRU;

*j*is the unitary calcium flux, i.e., the amount of Ca

_{Ca}^{2+}(in moles) that flows through a unit area per unit time if a single RyR is open;

*τ*is the time constant of calcium removal from the cytosol by the SERCA pump;

_{d}*D*stands for the effective diffusion coefficient of calcium ions that reflects Ca

_{Ca}^{2+}buffering by endogenous and exogenous buffers; and the vectors

**r**

*represent the positions of the sources. Spatial and temporal discretization of the sources is expressed via the Dirac delta function*

_{i}*δ*and the Heaviside step function

*H*, respectively. This means that the sources localized at

**r**

_{i}are activated at the moments

*t*. 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).

_{rd}The times *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 *n _{RyR}* channels at least one opens during Δ

*t*,

*P*, 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:

_{CRU}*P*is the probability that a single RyR channel opens during the interval Δ

_{RyR}*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

*k*represents the rate constant of the transition to the open state. Considering the relationship between

^{+}*k*, single-channel steady-state open probability

^{+}*P*(

_{O}*Ca*), and the average open time

*τ*, the probability that a single RyR opens within the interval Δ

_{o}*t*is

### Methods

#### Sensitivity of spontaneous spark properties to changes in cytosolic [Ca^{2+}].

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

^{2+}]. 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}*ϕ*/∂[Ca

_{spark}^{2+}] was introduced. The analytical expression for ∂

_{Ca}

*ϕ*was derived by differentiating Eq. 2, into which Eq. 1 has been substituted, with respect to [Ca

_{spark}^{2+}]. 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 ^{2+}]. Because calculation of *p _{O}* (Eq. S1), which is not defined by an analytical expression,

^{2+}] 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 *n _{RyR}* = 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 Ca

^{2+}with a flux

*J*for a period

_{Ca}*t*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

_{rd}*t*. The simulation was initialized by randomly assigning active, closed, and refractory state to individual CRUs with probabilities conforming to the values (

_{rp}*,*

^{ss}P_{spark}*,*

^{ss}P_{idle}*) given by Eqs. 3–6. The mean RyR open time was set to*

^{ss}P_{ref}*τ*= 15 ms (Tencerová et al., 2012; Gaburjakova and Gaburjakova, 2014), and the probability

_{o}*P*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 [Ca

_{O}^{2+}]. Afterward, the local dynamically calculated [Ca

^{2+}] at each CRU (see below) was used to calculate its

*P*(Eq. 10) during the current time step. The cytosolic concentration of Mg

_{CRU}^{2+}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,

*t*). The CRU parameters with their standard values are summarized in Table 2.

_{rp}CRUs in the active state are calcium sources with a flux *J _{Ca}*, and the released Ca

^{2+}ions diffuse throughout the domain with an apparent diffusion coefficient

*D*that accounts for calcium binding to endogenous and exogenous buffers in the cytosol. Reuptake of Ca

_{Ca}^{2+}ions from the cytosol into the SR by SERCA was represented by the time constant

*τ*. The flux

_{d}*J*is given by the number of open RyRs (

_{Ca}*n*) in the active CRU and the unitary flux

_{O}*j*through a single RyR (

_{Ca}*J*=

_{Ca}*n*×

_{O}*j*). The actual value of

_{Ca}*n*for each active CRU was set randomly according to the binomial probability distribution as explained in the Theory section (Eq. 7). The value of

_{O}*n*was generated using the “binornd” function from the Statistical Toolbox of MATLAB (version 8.3 [R2014a]; MathWorks) as follows:

_{O}*p*has been explained in the Theory section, and the procedure for

_{O}*p*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

_{O}*p*.

_{O}The calcium release flux through a single open RyR (the unitary flux *j _{Ca}*) was assigned by the following procedure. The spatiotemporal distribution during three-dimensional radial diffusion of Ca

^{2+}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 Ca

^{2+}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

*j*value was adjusted so that [Ca

_{Ca}^{2+}] at

*r*= 2 µm and

*t*= 15 ms was equal for both models for the standard value of the diffusion coefficient (

*D*= 30 µm

_{Ca}^{2}s

^{−1}). This procedure yielded the value

*j*= 0.481 µmol m

_{Ca}^{−2}s

^{−1}.

The time constant of Ca^{2+} sequestration by SERCA at low cytosolic [Ca^{2+}], *τ _{d}^{min}*, 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

*τ*is shown in Fig. S1. Importantly, the value of

_{d}*τ*is only very weakly cytosolic calcium dependent at [Ca

_{d}^{2+}] ≤ 0.25 µM, which justifies the use of Eq. 9 for calculation of the spatiotemporal evolution of calcium concentrations.

Model parameters that relate to Ca^{2+} 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 *t* = 3 ms. To assure a steady state for the whole analyzed time period, the simulation was started at a time *t* = −*t _{rp}* with a randomly assigned initial state of the array of CRUs as described above. At every time interval (

*t*,

_{i}*t*+ Δ

_{i}*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

*P*, calculated from the local value of cytosolic [Ca

_{CRU}^{2+}] at the given CRU at

*t*=

*t*(Eqs. 10 and 11). If the random number did not exceed the value of

_{i}*P*, the CRU became a Ca

_{CRU}^{2+}source with the flux

*n*×

_{O}*j*for the duration

_{Ca}*t*and subsequently became refractory for the duration

_{rd}*t*. The value of

_{rp}*n*was generated randomly after the binomial distribution (Eq. 12). The spatial distribution of cytosolic [Ca

_{O}^{2+}] in the domain at

*t*=

*t*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

_{i}*t*=

*t*− Δ

_{i}*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 [Ca

^{2+}] produced at a neighboring release site by calcium release flux of 0.98 pA, corresponding to

^{2+}] = 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,

### 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.

## RESULTS

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 Ca^{2+} 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 Ca^{2+} 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 (*D _{Ca}* = 30 µm

^{2}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

*n*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).

_{O}Visualization of a typical simulation under control conditions and at different cytosolic [Ca^{2+}] 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 [Ca

^{2+}] 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 [Ca

^{2+}] 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 [Ca

^{2+}], 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 [Ca

^{2+}] (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 *j _{Ca}* for an

*n*-fold increase of the parameter

*D*. Other parameters of the simulation remained unchanged. Visualization of a typical simulation at different values of

_{Ca}*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 [Ca

^{2+}], as apparent from the calcium dependence of the frequency of triggered calcium sparks (Δ

*ϕ*; 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

_{spark}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.

#### Theoretical considerations.

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 *P _{CRU}* = 0.003–0.015). To achieve high triggering probability at low spontaneous calcium spark frequency, the increase of local cytosolic [Ca

^{2+}] 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}

*ϕ*(see Materials and methods).

_{spark}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

Additionally, the calcium sensitivity of

#### Modulation of calcium release by changes in RyR gating parameters.

Calcium sparks were simulated in the range of cytosolic [Ca^{2+}] 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 Ca^{2+} and Mg^{2+} from the activation site (*K _{Ca}*,

*K*), the allosteric factors of Ca

_{Mg}^{2+}and Mg

^{2+}(

*f*,

_{Ca}*f*), and the probability of Mg

_{Mg}^{2+}dissociation from the activation site during a spark (

*p*)—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 [Ca

_{dis}^{2+}] 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

*K*,

_{Ca}*f*, and

_{Ca}*f*, an increase to 300% for

_{Mg}*K*, and an increase to 200% for

_{Mg}*p*. Fig. 3 (B and C) shows the values of the frequency of triggered sparks (Δ

_{dis}*ϕ*) and of the increase in the average number of open RyRs in triggered sparks

_{spark}^{2+}(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 Δ

*ϕ*and

_{spark}^{2+}] range, in contrast to data in Fig. 1. This result indicates that under control conditions, the tendency to form waves at elevated cytosolic Ca

^{2+}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 *K _{Ca}*,

*K*,

_{Mg}*f*, and

_{Ca}*p*, but not when

_{dis}*f*was decreased. Whereas modified

_{Mg}*K*,

_{Ca}*K*, and

_{Mg}*f*had the strongest effect on wave formation at low cytosolic [Ca

_{Ca}^{2+}] (Fig. 3 B) and resulted in an increased number of open RyRs in triggered sparks (Fig. 3 C), increased

*p*led to increased wave formation in the whole cytosolic [Ca

_{dis}^{2+}] 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 Δ

*ϕ*and

_{spark}_{Ca}

*ϕ*(Fig. 3 E, cf. panels corresponding to changed

_{spark}*K*,

_{Ca}*K*, and

_{Mg}*f*), and to the average number of open RyRs in spontaneous sparks

_{Ca}*p*). The calcium sensitivity of

_{dis}^{2+}] = 0.4 µM after a decrease of the parameters

*f*and

_{Ca}*K*, but not after an increase in

_{Ca}*K*(Fig. 3, cf. C and G).

_{Mg}Changes of the remaining parameters related to spontaneous spark formation—the closed state dissociation constant, *K _{O}*

_{00}(an inverse measure of RyR tendency to spontaneous opening), the Mg

^{2+}dissociation constant of the inhibition site,

*K*, and the refractory period

_{I}*t*—are examined in Fig. 4. Fig. 4 A shows the results of typical simulations for cytosolic [Ca

_{rp}^{2+}] 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

*K*

_{O}_{00}and

*t*and an increase to 300% for

_{rp}*K*. Neither of the changes brought about a strong tendency to wave formation (Fig. 4 A). The decrease of

_{I}*K*

_{O}_{00}and

*t*and the increase of

_{rp}*K*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).

_{I}### 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 (

*j*), and by the apparent calcium diffusion coefficient

_{Ca}*D*.

_{Ca}The effect of a threefold change of the parameters *D _{Ca}*,

*j*, and

_{Ca}*τ*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 [Ca

_{d}^{2+}] (Fig. 5 A, bottom row). In contrast to the data presented in Figs. 3 and 4, the highest values of Δ

*ϕ*for increased parameters were observed especially at intermediate [Ca

_{spark}^{2+}] (Fig. 5 B). The value of

*D*and

_{Ca}*j*, but above all for increased

_{Ca}*τ*.

_{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

*f*in Fig. 3 A and parameters

_{Mg}*K*

_{O}_{00},

*K*, and

_{I}*t*in Fig. 4 A) or by increased cytosolic [Ca

_{rp}^{2+}] (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}*D*,

_{Ca}*j*, and

_{Ca}*τ*), and the variables characterizing the triggered calcium release events (Δ

_{d}*ϕ*,

_{spark}The frequency of triggered sparks, Δ*ϕ _{spark}*, was dependent on the product

_{Ca}

*ϕ*and on the weighted average of the squares of variables characterizing calcium dynamics,

_{spark}*D*,

_{Ca}*j*, and

_{Ca}*τ*, resulting in the equation

_{d}*R*) = 0.92 and coefficient of variation (

^{2}*ĉ*) = 0.53 (solid lines in Figs. 3 B, 4 B, and 5 B). The proportionality of Δ

_{v}*ϕ*to

_{spark}_{Ca}

*ϕ*is caused by the commensurability between

_{spark}^{2+}] elevation at the site of the triggered spark. The lack of a direct dependence of Δ

*ϕ*on spontaneous spark frequency can be explained by the implicit proportionality between

_{spark}*ϕ*and ∂

_{spark}_{Ca}

*ϕ*. The quadratic dependence of Δ

_{spark}*ϕ*on the calcium-handling parameters (diffusion coefficient

_{spark}*D*, unitary calcium release flux

_{Ca}*j*, and time constant of calcium sequestration

_{Ca}*τ*) indicates nonlinear character of the spark-triggering process. The parameters A

_{d}_{1}–A

_{3}describing the relative effect of the variables

*D*,

_{Ca}*j*, and

_{Ca}*τ*are summarized in Table 4. The results are shown as solid lines in Figs. 3 B, 4 B, and 5 B.

_{d}The increase of the average number of open RyRs in triggered sparks, _{Ca}*ϕ _{spark}* and on a weighted function of the variables characterizing calcium dynamics,

*D*,

_{Ca}*j*, and

_{Ca}*τ*. The best description of the data was achieved by assuming a quadratic dependence of

_{d}*D*, and fourth-order dependence on the unitary calcium release flux

_{Ca}*j*and on the time constant of calcium sequestration

_{Ca}*τ*, resulting in the equation

_{d}*R*= 0.88 and

^{2}*ĉ*= 0.46 (solid lines in Figs. 3, 4, and 5 D). The higher order dependence of

_{v}_{0}describing the effect of the variable

_{1}–A

_{3}describing the relative effect of the variables

*D*,

_{Ca}*j*, and

_{Ca}*τ*are summarized in Table 5. The results are shown as solid lines in Figs. 3 C, 4 C, and 5 C.

_{d}Despite the simplicity of the equations and the large variety of the observed relationships between RyR gating parameters, cytosolic [Ca^{2+}], and the properties of triggered calcium sparks, the description of the calcium dependence of Δ*ϕ _{spark}* and

The relative effect of individual properties of spontaneous sparks (∂_{Ca}*ϕ _{spark}*,

*D*,

_{Ca}*j*, and

_{Ca}*τ*) 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 ∂

_{d}_{Ca}

*ϕ*and

_{spark}^{2+}. To suppress the effect of

*D*,

_{Ca}*j*, and

_{Ca}*τ*, the parameters of Eqs. 13 and 14—A

_{d}_{0}, A

_{1}, A

_{2}, and A

_{3}, 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 (Δ

*ϕ*) is controlled to a large extent by the calcium sensitivity of calcium spark frequency (∂

_{spark}_{Ca}

*ϕ*). The number of open RyRs in spontaneous sparks

_{spark}*p*, when

_{dis}*f*and

_{Ca}*K*in the low cytosolic [Ca

_{Mg}^{2+}] range. The relative effect of ∂

_{Ca}

*ϕ*and

_{spark}^{2+}]. At low cytosolic calcium, the effect of ∂

_{Ca}

*ϕ*may contribute up to 50% to

_{spark}^{2+}dissociation during the spark, which affects

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

*D*or

_{Ca}*j*is increased, its effect becomes similar in size or even larger than that of

_{Ca}*τ*. The relative contribution of the calcium-handling parameters

_{d}*D*,

_{Ca}*j*, and

_{Ca}*τ*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 (

_{d}*f*,

_{Ca}*p*), the contribution of

_{dis}*D*to the increase of the number of open RyRs in triggered sparks

_{Ca}*j*and

_{Ca}*τ*. Increase in any of the calcium-handling parameters strongly increased their effect on

_{d}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

*ϕ*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 Ca

_{spark}^{2+}; 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 Ca

^{2+}. For groups RyR-Act and CaH, there was a highly significant correlation (Pearson’s correlation coefficient [

*R*] = 0.90 and 0.89, respectively; P < 0.001) between Δ

_{P}*ϕ*and

_{spark}*ϕ*and

_{spark}*R*= 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

_{P}*t*, [Ca

_{rp}^{2+}] = 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 (Δ

*ϕ*= 0;

_{spark}*R*= 0.44; P = 0.18). Overall correlation between Δ

_{P}*ϕ*and

_{spark}*R*= 0.91; P < 0.001; slope, 4.88 ± 0.28 s). Thus, the parameter Δ

_{P}*ϕ*is suitable for quantification of the propensity to form calcium waves, and no additional organization of sparks into waves was observed in the simulations.

_{spark}The average number of sparks in the wave was relatively low *n _{Spark}* of individual waves (not depicted). However, even under conditions supporting the highest wave formation (decreased

*f*or increased

_{Ca}*τ*), 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 Ca

_{d}^{2+}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 (

*R*= −0.99; P < 0.002) and the fraction of waves terminated by refractoriness increased (

_{P}*R*= 0.99; P < 0.002) with cytosolic [Ca

_{P}^{2+}]. 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

_{Ca}

*ϕ*vs. by increased

_{spark}*K*

_{O}_{00}). The effect of increased

*K*

_{O}_{00}, 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

*K*), decreasing the magnesium affinity of the activation site (increase of the RyR gating parameter

_{Ca}*K*), increasing the positive allosteric effect of Ca

_{Mg}^{2+}binding to the RyR activation site (decrease of the RyR gating parameter

*f*), or increasing the probability of Mg

_{Ca}^{2+}unbinding (increase of the RyR gating parameter

*p*). 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

_{dis}*f*, and after combining the effect of decreased

_{Ca}*f*with the effect of increased

_{Ca}*K*

_{O}_{00}, is shown in Fig. 9 (A–C). The formation of waves arising because of a threefold decrease of

*f*(cf. Fig. 9, A and B) can be almost fully reversed by a threefold increase in

_{Ca}*K*

_{O}_{00}(Fig. 9 C). Statistical analysis for all tested RyR gating parameters is shown in Fig. 9 D.

The ability of changes in *K _{O}*

_{00}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.

## DISCUSSION

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 *K _{Ca}* or

*f*, or an increase in

_{Ca}*K*at 0.05 µM [Ca

_{Mg}^{2+}]; 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

*p*at 0.05 µM [Ca

_{dis}^{2+}]; Fig. 3); it may remain unchanged after an increase in the spontaneous spark frequency (e.g., a decrease in

*f*at 0.25 µM [Ca

_{Mg}^{2+}]; Fig. 3), but it may also decrease after an increase in the spontaneous spark frequency (e.g., an increase of [Ca

^{2+}] at

*f*= 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.

_{Ca}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 [Ca^{2+}] 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 (*K _{Ca}*,

*K*,

_{Mg}*f*). 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 (

_{Ca}*p*), a high propensity to form calcium waves could be observed in the whole examined range (0.05–0.4 µM) of cytosolic [Ca

_{dis}^{2+}].

### 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 [Ca^{2+}] 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 Ca^{2+} 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

_{Ca}

*ϕ*,

_{spark}*D*,

_{Ca}*j*, and

_{Ca}*τ*). 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 (∂

_{d}_{Ca}

*ϕ*) and of the number of open RyRs in spontaneous sparks

_{spark}#### RyR gating and wave termination.

Statistical analysis revealed that the frequency of triggered calcium sparks (Δ*ϕ _{spark}*) unequivocally determined the average length of waves

Wave termination was dependent on both cytosolic Ca^{2+} 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 Ca^{2+} 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.

### Limitations

### Gating model limitations

#### Luminal control.

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 *f _{Ca}*,

*f*, and

_{Mg}*K*

_{O}_{00}. Out of these parameters, only a change in the allosteric factor

*f*increases the propensity to form calcium waves.

_{Ca}#### Refractory period.

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 Mg^{2+} 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 Ca^{2+}, 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 Mg^{2+} under diastolic conditions and the finite rate of Mg^{2+} 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, *i _{Ca}*, was assumed constant throughout the spark. However, it can be assumed that

*i*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,

_{Ca}*i*stayed constant throughout the spark, whereas in the second simulation, it decreased linearly. The time integral of

_{Ca}*i*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

_{Ca}*i*does not affect the conclusions of this work.

_{Ca}#### One-dimensional approximation.

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.

#### Linear approximation.

The presented model does not involve binding of Ca^{2+} 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 *D _{Ca}* (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 [Ca^{2+}], as the cytosolic calcium-binding site of SERCA has submicromolar *K _{Ca}*. 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

*τ*(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 [Ca

_{d}^{2+}] (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.

### Physiological significance

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 Ca^{2+} binding to the RyR activation site and channel opening is increased, which can be described by a decrease of the parameter *f _{Ca}* in the aHTG model of RyR gating. Several arrhythmogenic RyR mutations (P2328S, Q4201R, V4653F) were shown to exhibit two to three times decreased Mg

^{2+}sensitivity (i.e., increased

*K*) 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.

_{Mg}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 [Ca^{2+}]-induced decrease of *f _{Ca}* 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 Mg^{2+}), 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 *i _{Ca}* 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 Ca

^{2+}(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 *K _{Ca}*,

*K*) and modulation of RyR activation by Ca

_{Mg}^{2+}in indirect fashion (changes in the allosteric coupling between Ca

^{2+}binding and RyR opening, defined by

*f*) 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 Mg

_{Ca}^{2+}dissociation constant (

*K*) versus its unbinding rate (

_{Mg}*p*) 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 (

_{dis}*K*,

_{Ca}*K*,

_{Mg}*f*, and

_{Ca}*p*). The results are reassuring, as inhibition of spontaneous RyR opening deeply and significantly suppressed the formation of calcium waves in all studied cases.

_{dis}### Conclusions

In brief, this study shows that changes in RyR gating may lead to the appearance of calcium waves at diastolic cytosolic [Ca^{2+}] 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 [Ca^{2+}] 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 (*K _{Ca}*,

*K*,

_{Mg}*f*,

_{Ca}*p*) and upon changes in cytosolic [Mg

_{dis}^{2+}]. 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}

*ϕ*, and parameters characterizing calcium dynamics of the cell (

_{spark}*D*,

_{Ca}*j*, and

_{Ca}*τ*).

_{d}## Acknowledgments

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.

## Footnotes

- Abbreviations used in this paper:
- aHTG model
- allosteric homotetrameric gating model
- CRU
- calcium release unit
- FDF
- fire–diffuse–fire
- SERCA
- 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/).