The dyads of cardiac myocytes contain ryanodine receptors (RYRs) that generate calcium sparks upon activation. To test how geometric factors of RYR distribution contribute to the formation of calcium sparks, which cannot be addressed experimentally, we performed in silico simulations on a large set of models of calcium release sites (CRSs). Our models covered the observed range of RYR number, density, and spatial arrangement. The calcium release function of CRSs was modeled by RYR openings, with an open probability dependent on concentrations of free Ca2+ and Mg2+ ions, in a rapidly buffered system, with a constant open RYR calcium current. We found that simulations of spontaneous sparks by repeatedly opening one of the RYRs in a CRS produced three different types of calcium release events (CREs) in any of the models. Transformation of simulated CREs into fluorescence signals yielded calcium sparks with characteristics close to the observed ones. CRE occurrence varied broadly with the spatial distribution of RYRs in the CRS but did not consistently correlate with RYR number, surface density, or calcium current. However, it correlated with RYR coupling strength, defined as the weighted product of RYR vicinity and calcium current, so that CRE characteristics of all models followed the same state-response function. This finding revealed the synergy between structure and function of CRSs in shaping dyad function. Lastly, rearrangements of RYRs simulating hypothetical experiments on splitting and compaction of a dyad revealed an increased propensity to generate spontaneous sparks and an overall increase in calcium release in smaller and more compact dyads, thus underlying the importance and physiological role of RYR arrangement in cardiac myocytes.

Calcium spark is a miniature transient increase of fluorescence intensity observable in the cytosol of cardiac myocytes using calcium-sensitive fluorescent dyes (Cheng et al., 1993; López-López et al., 1994). It represents a local event of calcium release that appears at dyads under specific experimental conditions. The dyads are intracellular structural complexes of cardiac myocytes, where terminal cisternae of the SR connect to the sarcolemma (Sun et al., 1995; Franzini-Armstrong et al., 1998). Under physiological conditions, the local calcium release events (CREs) are triggered in synchrony at all dyads in the whole volume of a myocyte in response to the influx of Ca2+ ions through sarcolemmal calcium channels activated by incoming action potential (Cannell et al., 1995; López-López et al., 1995). The resulting global calcium release provides myofibrillar sarcomeres with calcium ions needed for contraction (Eisner et al., 2017).

Calcium sparks also occur in resting cardiac myocytes. Under physiological conditions, these spontaneous sparks occur at low frequency; however, under pathological conditions, their frequency may increase and even initiate arrhythmogenic calcium waves (Cheng et al., 1996; Lukyanenko and Gyorke, 1999). Spontaneous calcium sparks also originate at dyads, but independently of the sarcolemmal voltage-activated calcium channels (Cheng et al., 1993). Their activation results from a spontaneous opening of one of the ryanodine receptors (RYRs), calcium release channels clustered at a terminal cisterna. The released Ca2+ ions activate nearby RYRs by the calcium-induced calcium release mechanism (Fabiato, 1983) and thus give rise to the local spark. Calcium sparks depend primarily on RYR function, which is regulated by cytosolic and luminal Ca2+ ions, cytosolic Mg2+, and cytosolic ATP2−, as well as by numerous regulatory proteins and processes (Meissner, 2002, 2004). For physical reasons, it should also depend on the group characteristics of RYR clusters. The question of how the distribution of RYRs in the cisternae affects the properties of calcium sparks remains open.

The structure of terminal cisternae studied by EM revealed common features in their construction (Franzini-Armstrong et al., 1999; Hayashi et al., 2009) as well as their substantial structural variability (Novotova et al., 2009; Novotova et al., 2013; Novotová et al., 2020). Super-resolution localization studies revealed a tendency of RYRs to form irregular clusters of variable size and surface density (Baddeley et al., 2009; Hou et al., 2015; Macquaide et al., 2015; Galice et al., 2018; Jayasinghe et al., 2018; Kolstad et al., 2018; Asghari et al., 2020). Changes in the size, shape, and content of terminal cisternae, regarding the number of RYRs and amount of calsequestrin, accompany changes in the structure of dyads (Asghari et al., 2020; Drum et al., 2020; Novotová et al., 2020). It should be noted that the positions of spontaneous calcium sparks, termed calcium release sites (CRSs), were associated with terminal cisternae according to indirect observations, such as colocalization of the spark fluorescence with that of t-tubules and z-discs (Parker et al., 1996; Cleemann et al., 1998; Hiess et al., 2015). Our understanding of the terminal cisterna as the substrate of calcium sparks is complicated by the difficulty to estimate explicitly the composition of the source of calcium flux. It is not clear how RYRs of a terminal cisterna participate in the formation of calcium sparks. Is it the whole cisterna, one of its RYR clusters, or even only a single RYR, which determines calcium sparks?

The calcium release function of CRSs can be studied only indirectly, in particular by monitoring CREs, such as calcium sparks or calcium spikes (Cheng et al., 1993; Song et al., 1998). Despite the recent progress in technologies, these approaches have limited spatiotemporal resolution and signal-to-noise ratio, inherent to the calcium-sensitive fluorescent dyes (Smith et al., 1998; Zahradníková et al., 2007). CREs vary broadly from invisible events, named quarks (Lipp and Niggli, 1996), through miniature “quarky calcium release” (Brochet et al., 2011; Shang et al., 2014), to brief sparks (Cheng et al., 1993) and prolonged sparks (Zima et al., 2008; Brochet et al., 2011). The CRE signal varies even if recorded repeatedly at the same spot (Sobie et al., 2006; Zima et al., 2008; Janíček et al., 2012; Shang et al., 2014). Interpretation of local fluorescent signals is difficult since it is based on data obtained either indirectly or in artificial systems. Therefore, theoretical approaches could be of help to formulate and verify working hypotheses (Soeller and Cannell, 2004; Williams et al., 2007; Soeller et al., 2009; Williams et al., 2011; Stern et al., 2013).

According to recent understanding, calcium sparks result from group activity of numerous RYRs clustered at CRSs. A close association of RYRs was observed in cardiac dyads (Asghari et al., 2014) and purified cardiac RYRs (Cabra et al., 2016); however, only a small fraction of RYRs, whether in situ or in vitro, physically interact (Zahradníková et al., 2010a; Jayasinghe et al., 2018). The first estimates of the junctional SR surface area per RYR were close to 900 nm2 (i.e., 30 × 30 nm; Franzini-Armstrong et al., 1999; Soeller et al., 2007). This corresponds to a tightly packed RYR cluster. The more recent estimates report values in the range of 3,000−12,500 nm2 (Baddeley et al., 2009; Hayashi et al., 2009; Hou et al., 2015; Galice et al., 2018; Jayasinghe et al., 2018). The average number of RYRs per CRS, often used as a measure of CRS size, was estimated in the range of 8–22 for terminal cisternae attached to the surface sarcolemma (Baddeley et al., 2009; Jayasinghe et al., 2018), and 8–100 for terminal cisternae attached to the tubular sarcolemma (Soeller et al., 2007; Hayashi et al., 2009; Scriven et al., 2010; Hou et al., 2015; Galice et al., 2018). The existence of clusters with RYRs separated by the edge-to-edge gap of >100–150 nm led to the postulation of superclusters (Baddeley et al., 2009; Hou et al., 2015; Macquaide et al., 2015). Thus, the definition of CRS size by the number of its RYRs is somewhat fuzzy. However, the question of how the CRS size, RYR number, and RYR surface density are related to calcium spark amplitude and duration is of fundamental physiological importance, for instance, for understanding the changing calcium release function during postnatal development (Ziman et al., 2010; Louch et al., 2015; Macková et al., 2017) or in various pathologies (Song et al., 2006; Wu et al., 2012).

The understanding of cardiac dyads can be advanced by mathematical modeling, which allows for the exploration of the whole range of experimental data through the parametric space of defined models. In this modeling study, we addressed mainly the issues of calcium spark variability related to geometrical factors of RYR arrangement at CRSs. We developed a set of models to compare CREs produced by CRSs differing in RYR number, arrangement, surface density, and single-channel calcium current. The model parameters spanned the whole range of relevant findings published recently. Activation of RYRs was simulated by a gating model conforming to the experimentally observed sensitivity of the cardiac RYR to activation by Ca2+ ions and inhibition by Mg2+ ions (Laver et al., 1997; Zahradníková et al., 2003). This allowed us to elucidate how the CRS determinants affect the properties of local CREs.

### Simulations

Programs for construction of the topology of CRSs and simulation of the activity of RYRs at CRSs were written in C++. The algorithms are described in Model formulation. Both programs were parallelized using the OpenMP API (https://www.openmp.org). Parallel generation of random numbers was performed using OpenMP Parallel Statistical Random Number Generator OMPRNG (M. Bognar, University of Iowa, Iowa City, IA; https://homepage.stat.uiowa.edu/~mbognar/omprng/). The Uran supercomputer at Krasovskii Institute of Mathematics and Mechanics of the Ural Branch of the Russian Academy of Sciences was used for part of the model calculations. Cytosolic calcium signals were simulated using CalC (Matveev et al., 2002).

### Analyses

Data analysis was performed in OriginPRO version 2020b (OriginLab). The data were fitted by theoretical functions using the Levenberg-Marquardt method.

### Model formulation

Models of a CRS were represented by a group of RYRs characterized by geometrical and functional descriptors. The geometrical descriptors included the number of RYRs in the CRS (NRYR), the area of CRS per one RYR (AC/R), arrangement patterns of RYR distribution (checkerboard [ChB], side by side [SbS], and isolated RYRs), and coordinates of RYRs. The functional descriptors included the RYR single-channel calcium current (iCa), the on- and off-rate constants of RYR gating (kco and koc, respectively), and parameters of calcium diffusion and binding to buffers. The values of NRYR, AC/R, and iCa, as well as arrangement patterns were varied within the range of published values, and their impact on the generation of CREs was evaluated.

### The model of RYR gating

We used a simple two-state RYR gating model to minimize requirements for computational power (Eq. 1), with the open probability dependent on free Ca2+ and Mg2+ concentration:
$C⇌kockcoO.$
(1)
The gating model obeyed experimental observations—namely, the RYR open probability at low Ca2+ concentration and the effect of Mg2+ binding to the RYR inhibition sites (Laver et al., 1997; Zahradníková et al., 2003; Tencerová et al., 2012). The value of the off-rate constant koc = 0.45 ms−1 was set in correspondence to the mean RYR open time (2.2 ms), estimated in the presence of ATP and Mg2+ on the cytosolic RYR side and at a realistic Ca2+ concentration on the luminal RYR side (Zahradníková et al., 2003). Similar values under similar conditions were reported by Györke and Györke (1998) and Laver and Honen (2008). The on-rate constant kco was made dependent on the cytosolic [Ca2+] and [Mg2+] as
$kco([Ca2+],[Mg2+])=kocPa([Ca2+],[Mg2+])1−Pa([Ca2+],[Mg2+])$
(2)
(Iaparov et al., 2019), where Pa is the Ca2+- and Mg2+-dependent open probability of the allosteric homotetrameric gating model according to Zahradníková et al. (2010b). The [Ca2+] and [Mg2+] are free ion concentrations calculated for coordinates of the center of an RYR. The steady-state open probability of this model,
$PO=kcokco+koc,$
(3)
was thus dependent on [Ca2+] and [Mg2+], as shown in Fig. 1 for 1 mM [Mg2+] used in all CRS simulations. Under these conditions, the maximum steady-state open probability $POmax$ was 0.31.

The variable clustering of RYRs in terminal cisternae was modeled by generating a large set of CRSs with different RYR arrangement that was varied following experimental data. Four topologically different CRS arrangement patterns (Fig. 2 and Table 1) were made by placing 10, 20, 40, or 80 RYRs on the CRS plane at calculated or experimentally determined positions (see below) to encompass the range of reported CRS sizes and RYR surface densities (Baddeley et al., 2009; Hou et al., 2015; Galice et al., 2018; Jayasinghe et al., 2018).

To characterize the packing of RYRs in CRSs, we used the AC/R (nm2), that is, the average area of CRS associated with one RYR, which is more suitable for geometrical considerations than the RYR surface density (nm−2). The NRYR and the AC/R were used to calculate the side a of the CRS as $a=NRyR.AC/R.$ Three CRS density type models were used in this study to cover the whole range of reported AC/R values.

The compact-type CRSs were of two variants. The square variants of compact CRSs (Fig. 2 A) contained the given number of RYRs placed side by side on a rectangular grid of 30 × 30-nm squares. To accommodate the given number of RYRs in the square grid, the empty positions were placed at random. The resulting AC/R was 1,440, 1,125, 1,103, and 911 nm2 per RYR, respectively, for 10, 20, 40, and 80 RYRs in the CRS (areas equivalent to squares with sides of 38, 34, 33, and 30 nm). For each NRYR, a set of five compact square CRSs was generated (altogether 20 CRSs). The irregular variants of compact CRS (Fig. 2 B) used the experimentally determined coordinates (Walker et al., 2015) of 9–96 RYRs placed on a contact network grid with 31-nm spacing. Altogether, 107 of these CRSs had AC/R equal to 961 nm2 per RYR (the area equivalent to a square with a side of 31 nm).

The narrow-type CRSs (Fig. 2 C) had RYR surface density corresponding to that estimated by Jayasinghe et al., 2018. All narrow CRSs had an AC/R of 3,481 nm2 (the area equivalent to a square with a side of 59 nm). For each NRYR and each pattern (Table 1), a set of 100 CRSs with different RYR coordinates was generated (altogether 1,200 narrow CRSs).

The wide-type CRSs (Fig. 2 D) had RYR surface density corresponding to that estimated by Galice et al., 2018. All wide CRSs had an AC/R of 12,544 nm2 (the area equivalent to a square with a side of 112 nm). For each NRYR and each pattern (Table 1), a set of 100 CRSs of different RYR coordinates was generated (altogether 1,200 wide CRSs).

Arrangement of RYRs in compact CRSs followed the side-by-side configuration (Yin et al., 2005; Groff and Smith, 2008; Walker et al., 2015). The topology of the narrow and wide CRSs consisted of RYRs arranged into three patterns with defined proportions of neighbors in the ChB configuration, the SbS configuration, and the isolated position, as observed by EM (Asghari et al., 2014) under different experimental conditions (Table 1). Pattern A was dominated by RYRs in the ChB configuration, typical for low cytosolic Mg2+ concentrations. Pattern C was dominated by RYRs in the SbS configuration, typical for high cytosolic Mg2+ concentrations. The balanced pattern (B) contained the ChB, SbS, and isolated RYR configurations, as observed for near-physiological conditions. The placement of RYRs in the narrow and wide CRSs according to patterns A, B, or C was solved as an optimization problem (Iaparov et al., 2017; Iaparov et al., 2019). The prescribed distribution was achieved with the use of the global optimization algorithm GSA (genetic simulated annealing; Xu et al., 2011), which is a combination of a genetic algorithm and simulated annealing. Briefly, the algorithm generates a random distribution of the given number of RYRs (NRYR) and then minimizes the sum of squares of differences between the prescribed and the achieved number of channels in each arrangement. Minimization proceeded until the prescribed configuration was exactly achieved. The GSA is a stochastic algorithm; thus, each execution generated a CRS with different RYR coordinates but with the same distribution of the configurations (Table 1). The RYR coordinates of individual CRSs were used to calculate their set of inter-RYR distances.

### Simulation and analysis of CREs

Each CRS model was subjected to simulations of CREs in which the stochastically simulated activity of all RYRs was recorded for 100 ms. The stochastic nature of RYR gating results in temporal fluctuations of open RYRs; therefore, each simulation of the CRS activity was repeated 10,000 times for each CRS model. This means that each RYR of the CRS was initialized by 10,000/NRYR times (on average).

To unveil the impact of calcium current passed by a single RYR on spark generation by CRS models of diverse types, we used four iCa values of 0.04, 0.15, 0.25, and 0.4 pA. These were kept constant during simulations. The calcium current through a single RYR or a single CRS has not yet been determined in situ. Numerous experimental difficulties make the estimates ambiguous because of their dependence on many variables that are difficult to control in experiments. Janíček et al. (2012) concluded that the calcium release flux per single RYR should be near 0.4 pA, a value satisfying estimates by independent experiments. In previous CRS modeling studies, iCa of 0.04 pA (Groff and Smith, 2008), 0.15 pA (Walker et al., 2015), or 0.6 pA (Cannell et al., 2013) were used for spark simulations.

The concentration of calcium ions inside the CRS is not homogenous (Valent et al., 2007; Walker et al., 2014). At each locality, it depends on the positions of open RYRs and the single-channel calcium current in a complex manner, with it being affected by the diffusion rates, the geometry of the CRS, buffering by mobile and immobile buffers, and calcium-binding kinetics. In this study, calcium concentrations at each RYR position were calculated using the Naraghi and Neher approximation for calcium diffusion in the cytosolic half-space delimited by the plane of the CRS, based on the supposition of instantaneous formation of a steady-state calcium gradient upon RYR opening and sufficient calcium buffering capacity of the cytosol (Naraghi and Neher, 1997). The spatial [Ca2+] profile around an RYR is then an explicit function of mean reaction times of calcium ions with buffers and of the buffering capacity and mobility of each buffer. Buffering of Ca2+ by calmodulin, ATP, and Fluo-4 was calculated using parameters of Ca2+ diffusion and the concentrations, equilibrium constants, and rate constants of Ca2+ binding to ATP, calmodulin, and Fluo-4 taken from Hake et al. (2012). A record of the number of open RYRs in time was generated using the Gillespie algorithm (Gillespie, 1977). Briefly, each CRE was started by an opening of the selected RYR. Calcium influx changed calcium distribution in the CRS and, correspondingly, the opening rate constants of individual RYRs. These stayed constant until the next opening or closure of any RYR in the CRS because of the instantaneous and stationary calcium gradient approximation (see above). The Gillespie algorithm determines the channel that changes its state—opens or closes—and the time of this change. The identity index of the RYR that would change its state was sampled from the distribution of probabilities of state transitions of individual RYRs, which was for each RYR equal to the ratio of its transition rate constant to the sum of rate constants over all RYRs. The time of the change was sampled from the exponential distribution of waiting times calculated for a rate constant equal to the sum of instantaneous values of transition rate constants over all RYRs in the CRS (kco for closed RYRs and koc for open RYRs). The resulting transition led to a new calcium distribution in the CRS and a new set of opening rate constants. The algorithm was repeated until the end of the simulation.

The simulated traces of RYR activity (i.e., the CREs), were characterized by the number of simultaneously open RYRs at the peak of the trace $NOpeak,$ by the time to peak (TTP; the time to the first occurrence of $NOpeak),$ and by the relative amplitude of the peak (Arel), normalized to the number of RYRs in the CRS as:
$Arel=(NOpeak−1)/NRyR,$
(4)
where $NOpeak-1$ is the maximal number of simultaneously recruited RYRs. These CRE characteristics were used for constructing histograms and for further analyses as described below.

### Calculations of fluorescence signals from CRE traces

The fluorescence signal was calculated from the time course of the total calcium current produced by a CRE taken as a spherical source with a 200-nm diameter and considering spherically symmetrical calcium diffusion, binding of calcium to cytosolic buffers, and fluorescence proportional to the concentration of calcium-bound Fluo-4, as described in Hake et al. (2012) and Walker et al. (2014). Parameters of binding and diffusion of ATP, calmodulin, troponin, and Fluo-4 were taken from Hake et al. (2012), and those of the sarcolemmal and SR membrane buffers were taken from Zahradníková et al. (2007). The calculated spatial profile of Ca2+-bound Fluo-4 concentration was convolved with a Gaussian point spread function to approximate the image observed by fluorescent microscopy (Smith et al., 1998).

According to the calcium-induced calcium release theory, Ca2+ ions are the interaction element between RYRs that spreads by diffusion. In published models of CRSs, this fact was handled in different ways. For CRSs with uniform RYR arrangement, Williams et al. (2011) introduced the mean-field approximation of Ca2+-mediated RYR–RYR interaction. This approximation does not apply to the general case of a nonuniform arrangement of RYRs (Iaparov et al., 2019). For compactly packed RYRs on a regular lattice, Walker et al. (2014) used the maximum eigenvalue of the RYR adjacency matrix as an overall parameter characterizing RYR arrangement. The adjacency principle, however, does not adequately characterize the geometry of nonadjacent freely distributed RYRs. Other measures of RYR arrangement, namely the average nearest-neighbor distance, the four nearest-neighbor distances, and the mean RYR occupancy, all related to the RYR surface density and mean RYR–RYR distance in CRSs (Cosi et al., 2019), correlated with parameters of simulated calcium sparks relatively weakly. Recently, we have found (Iaparov et al., 2017) that a promising predictor of spark occurrence for general CRSs could be the maximum eigenvalue of the RYR inverse distance matrix, which, however, is difficult to perceive for its somewhat abstract meaning. Therefore, we define here new geometric parameters, namely, the closeness and vicinity that characterize the spatial disposition of RYRs in a general CRS more conceivably and allow for interpretation of simulated CREs on common grounds, as demonstrated in Results.

Let us consider the steady-state solution of the diffusion equation that relates the Ca2+ distribution to the positions of open RYRs. Generally, in an unbuffered system, the concentration of Ca2+ ions at individual RYR positions is proportional to the sum of reciprocal values of distances Dij to the open RYRs, [Ca2+]RYRi ∝ ∑j1/Dij (Crank, 1975; nota bene: the presence of calcium buffers, considered in this study, decreases the concentration of free Ca2+ ions and distorts their steady-state distribution [Naraghi and Neher, 1997]; thus, the above relationship between RYR positions and Ca2+ distribution should be understood as an upper bound in formulations of Ca2+ mediated RYR–RYR interactions used below). The inverse of the distance can be thought of as the closeness Cij = 1/Dij, where indices i and j correspond to elements of the RYR distance matrix. This means that the formation of a CRE by interacting RYRs may be related to the RYR closeness matrix C composed of elements Cij, where the elements of the ith row correspond to the closeness of the remaining RYRs (ji) to the ith RYR. Let us further consider a linear transformation on a CRS plane yielding all coordinates of its elements multiplied by a factor K. This transforms the distances between RYRs in proportion to K, and the surface area per RYR in proportion to the square of K. Thus, when the surface area of a CRS, given as SCRS = AC/R × NRYR, is changed, the distances between RYRs are changed proportionally to the square root of the change in surface area—that is, Cij$1/AC/R$ when NRYR is constant, or Cij$1/NRyR$ when AC/R is constant. Considering these relationships allows us to characterize RYR distribution in the CRS by a new parameter—the normalized closeness $vij=CijNRyR$—which we named the vicinity of ith RYR to jth RYR. The effect of Ca2+ flux through the ith RYR on the whole CRS can then be considered a function of the average vicinity vi of the ith RYR to the NRYR − 1 remaining RYRs defined as
$vi=∑j,j≠ivijNRyR−1=NRyR∑j,j≠iCijNRyR−1.$
(5)
The average vicinity, in contrast to the average closeness, does not change when the surface area of CRS is altered by changing NRYR at a given AC/R. The whole CRS can be then characterized by the mean vicinity of all its RYRs (i.e., the group RYR vicinity v):
$v=∑iviNRyR.$
(6)
The group RYR vicinity defined in this way thus represents a CRS state variable characterizing the spatial component of the RYR–RYR interaction potential. As such, it will be used later for the analysis of simulated CREs. In contrast to the surface density, the vicinity characterizes the distribution of RYRs in the CRS.
Since CRS activity upon the opening of the ith RYR depends also on the calcium current through this RYR, we defined a combined state variable φi pertinent to and characterizing the ith RYR as the weighed product of calcium current and RYR vicinity:
(7)
where the power α is the weight factor for relative contributions of calcium current and vicinity of the ith initiating RYR to the CRS response. Developing this idea to account for CRS activation by the opening of any of its RYRs, the group RYR vicinity v (Eq. 6) allows us to define a state variable for the whole CRS as
(8)
From the point of physical meaning, the CRS state variables φi and φ could be considered to be indicators of RYR coupling strength, since they increase with increasing vicinity and with increasing calcium influx.
Values of the exponent α were found by approximating the relations of CRE characteristics to their respective RYR coupling strengths—φi or φ, as appropriate—using a general state-response function:
(9)
where x stays for φi or φ, Fmax is the maximal response, K is x at half-maximum, and n is the slope factor.

### Model limitations

Several simplifications were applied to keep computational costs reasonable. The model of RYR gating provides a realistic calcium dependence of RYR open probability; however, both the allosteric effect of Ca2+ and Mg2+ binding on RYR activation (Zahradník et al., 2005; Zahradníková et al., 2010b) and the Mg2+ binding at the inhibition site (Laver et al., 1997; Zahradníková et al., 2003) are modeled by a decrease of the RYR activation rate. As a result, the kinetics of RYR activation may differ to some extent.

The models of CRSs do not account for the effect of SR lumen depletion on single-channel calcium current. This simplification was necessary to understand the effect of geometric factors alone. Nevertheless, the effect of ica amplitude on the formation of CREs was tested extensively.

The use of a simplified model of buffered calcium diffusion in semi-infinite space (Naraghi and Neher, 1997), instead of a full reaction-diffusion model of the dyadic gap (Valent et al., 2007), substantially speeded up calculations but also undervalued the buildup of calcium concentrations at individual RYRs. However, considering that the activation of RYRs after an RYR opening proceeds much slower than the build-up of the calcium gradient (Soeller and Cannell, 1997; Valent et al., 2007), this simplification should not weaken the correlations described in Results.

### Online supplemental material

The activity of individual RYRs and the corresponding calcium concentration changes in the CRS plane, corresponding to the records of RYR activity in Fig. 3 A, are shown in the videos. Video 1 shows the temporal evolution of RYR activity and calcium concentration changes during a typical quark. Video 2 shows the temporal evolution of RYR activity and calcium concentration changes during a typical blip. Video 3 shows the temporal evolution of RYR activity and calcium concentration changes during a typical spark.

This study surveys CREs generated by 10,108 CRS models. Each model CRS was defined by a set of CRS determinants covering the range of published data, including the number of RYRs, CRS density type, the arrangement of RYRs into specific patterns, and single-channel calcium current values. CRE simulation was initiated by opening of one of the RYRs in the CRS and recording for 100 ms. CRE records were characterized by their $NOpeak,$ relative amplitude, TTP, and the molar amount of calcium ions released during the event (nCa). CRS activity was characterized by the average relative amplitude Ārel and the fraction of different response types and their average time course (see below) calculated from 10,000 simulations.

### Basic characteristics of simulated CREs

CRSs responded to the initial RYR opening by three very distinct types of CREs (Fig. 3), classified and named as follows.

The smallest CRE responses consisted of the initial RYR opening alone, which did not trigger other RYR openings. We named them quarks in analogy to the smallest CREs experimentally observed in situ (Lipp and Niggli, 1996; Brochet et al., 2011; Wescott et al., 2016). Durations of simulated quarks followed exponential distribution. The fluorescence image of a quark (Fig. 3 A) showed peak intensity below experimental resolution.

Intermediate CRE responses, emerging as openings of a few RYRs that deactivated rapidly, were named blips. At the peak of blips, the number of simultaneously open RYRs was typically between two and eight. Typical blips peaked in <1.5 ms and ceased completely in <15 ms. The fluorescence image of a blip (Fig. 3 A) showed peak intensity near the limit of experimental resolution (nota bene: the term blip used here is not meant to refer to calcium transients described by Swillens et al. [1998] as blips for inositol 1,4,5-trisphosphate receptor activity).

Responses containing openings of many RYRs with TTP typically above 3 ms that deactivated typically within 50 ms were named as sparks for their obvious similarity to real calcium sparks. The fluorescence image of a spark (Fig. 3 A) showed peak intensity above the limit of experimental resolution at all inspected iCa values.

The simulated gating activity of individual RYRs and the corresponding changes of calcium concentration in the CRS plane for quarks, blips, and sparks are demonstrated in Video 1, Video 2, and Video 3.

The amplitude histograms of all CREs (see the example in Fig. 3 B) showed a clear nadir under most CRS settings. The histograms were used to estimate the relative occurrences (fractions) of quarks, blips, and sparks in the simulation data sets obtained for all CRS determinants. The fraction of CREs in the first bin $NOpeak=1$ was assigned to quarks (Fq). The sum of CRE fractions in bins from $NOpeak=2$ to the nadir of the histogram was assigned to blips (Fb), and the sum of CRE fractions for bins above the nadir was assigned to sparks (Fs). The TTP histograms of simulated CRSs (Fig. 3 C) usually showed a monotonic decrease.

Correlations of peak amplitudes and TTPs showed that quarks, blips, and sparks could be distinguished from each other as regions of distinct fractional occurrence, as exemplified for selected CRS models in Fig. 4. The quarks are unique for their constant amplitude and TTP $(NOpeak=1;$TTP = 0) by definition. The blips have a TTP almost independent of their amplitude, while the TTP of sparks varies the broader, the higher is their amplitude. Similarly, the TTP of blips was shorter in CRSs of higher iCa (top row), more compact density type (middle row), and higher NRYR (bottom row). In the case of sparks, the $NOpeak$ and the fractional occurrence were higher and the TTP was longer in CRSs of higher iCa, more compact density type, and higher NRYR. This apparent tripartite CRS activity stems solely from the stochastic gating of RYRs. Its relation to the CRS determinants will be analyzed below.

### CREs depend on CRS determinants

The distributions of $NOpeak$ values varied substantially with CRS determinants (Fig. 5). The quarks and blips prevailed when the RYR number, density, and iCa were low, while the sparks dominated when at least one of these parameters was high. When iCa was very low, the sparks activated exclusively in large compact CRSs, but rarely (not shown). In CRSs with NRYR = 10, the distribution of $NOpeak$ values often did not have a clear nadir; thus, the distinction between blips and sparks was not so prominent.

The $NOpeak$ in individual sparks never reached NRYR, but due to the occurrence of approximately four reopenings per RYR on average, it often highly exceeded the expected values given by NRYR × $POpeak.$$NOpeak$ approached NRYR only in the smallest compact CRSs. In the largest compact CRSs, $NOpeak$ reached only ∼55% of NRYR, even at the largest iCa values. Generally, lower iCa values or lower vicinities resulted in a lower fraction of RYRs open at the peak of the CRE. This was due to the lower [Ca2+] at RYRs distant from a source RYR, resulting in their lower probability of opening and longer closures (Fig. 2, and Videos 1, 2, and 3).

Ca2+-mediated interactions between RYRs depend on distances between RYRs. The distances have a different distribution in different CRSs of the same RYR surface density. The distribution of RYRs in CRSs is conveniently characterized by their group vicinity (Eq. 6), which increases with the extent of RYR clustering. CRS models of the three density types had group RYR vicinities different to the extent that their ranges almost did not overlap (compact CRSs: v = 0.055 − 0.87; narrow CRSs: v = 0.031 − 0.051; wide CRSs: v = 0.021 − 0.036; see Table 2 in Iaparov et al. [2020] Preprint for details).

CRSs of identical density and pattern type, but of different RYR distribution, yielded different amplitude distributions of their simulated CREs (Fig. 6). The probability that the two CRE distributions are the same for the two CRSs differing only in their RYR group vicinity v was extremely low. Thus, the different fractional occurrence of quarks, blips, and sparks in the two CRSs was highly statistically significant. This result can be seen as evidence that the distribution of RYRs in real dyads affects calcium release activity. This also means, however, that the RYR group vicinity is a good integral descriptor of RYR clustering.

### CRS activity relates to CRS coupling strength quantitatively

The above analysis confirmed a close relation of CRS activity to CRS determinants (Figs. 4, 5, and 6); however, the characteristics of CREs did not correlate with CRS determinants, as exemplified by the average relative amplitude of CREs in Fig. 7 A. A similar lack of correlation was present in the relative occurrence of different CRE types. The only exception was a partial correlation of CRS activity with RYR group vicinity v, albeit different at different iCa values (Fig. 7 A, bottom). Therefore, we correlated the CRE characteristics with the RYR coupling strength φ that combines the RYR group vicinity and iCa (Eq. 8). The weight exponent α in Eq. 8 was found by approximating the relation of estimated CREs characteristics (Arel, Fq, Fb, Fs) to the coupling strength φ of their respective CRS models, defined by Eq. 8 using state-response functions (Eq. 9). Since Arel and Fs increased with increasing φ, both could be approximated by Eq. 9 directly. However, as the relationships between Fq or Fb and φ were decreasing or bell-shaped functions, we applied Eq. 9 through complementary relationships based on the fraction of non-quarks Fnq(φ) = 1 − Fq(φ); that is, Fq(φ) = 1 − Fnq(φ) and Fb(φ) = Fnq(φ) − Fs(φ). The whole dataset of the evaluated CRE characteristics, paired with their corresponding φ values (Eq. 8), was fitted collectively using a common value of the weight factor α. This was possible because α is related to the cytosolic calcium buffering, which is homogeneous within the CRS.

The CRE characteristics correlated very well with RYR coupling strength φ for α = 0.67 (Fig. 7 B). Agreement between the data and the theoretical functions was very high (Table 2). Of note is that the values of φ varied 25 times across all wide, narrow, compact, and contact network CRS models.

The existence of a functional relationship between φ and CRE characteristics for all CRS density types and the whole range of NRYR and iCa values proves that φ is a true descriptor of the CRS state. Thus, the coupling strength φ predicts the propensity of a CRS to produce sparks upon a spontaneous opening of one of its RYRs.

### Kinetics of CREs and the amount of calcium released relate to CRS determinants

Hand in hand with the amplitudes and relative occurrences, the kinetics of CREs also varied. For convenience, we analyzed the averaged CREs since individual CREs fluctuated substantially. A single quark terminates abruptly upon closure of the open RYR at a random time; however, the averaged quarks declined exponentially (Fig. 8 A, top). The rate constant of the average quark decline was faster than the off-rate constant of the RYR gating model since RYR openings in CRSs were curtailed by activation of blips and sparks. As a result, a typical decay time of quarks was ∼2.2 ms for small, wide CRS and ∼0.2 ms for large, compact CRSs. This also means that larger and denser RYR clusters lead to faster responses of blips and sparks to the initial RYR opening.

The time course of averaged blips was always biphasic (Fig. 8 A, middle). Their TTP was ∼1–1.5 ms in CRSs of low φ values, but at large φ values, it dropped to 0.1 ms due to the increased rate of spark formation (not shown). Typically, the blips terminated in full within 5–20 ms.

The time course of averaged sparks depended strongly on CRS determinants (Fig. 8 A, bottom). For the narrow CRS and iCa of 0.25 pA, the rate of spark activation was remarkably slower when the number of RYRs was larger. For 10 and 20 RYRs in the CRS, the averaged sparks peaked between 3 and 7 ms, while for 40 and 80 RYRs, the averaged sparks peaked between 10 and 30 ms. This indicates that even when the formation of the Ca gradient is instantaneous, as defined by the diffusion model, the spread of RYR activation in the CRS takes its time (see Videos 1, 2, and 3).

Termination of sparks was very complex. During individual sparks, the fraction of open RYRs fluctuated strongly and sometimes highly exceeded the steady-state PO. In large CRSs of high φ values, the average spark did not decline markedly until the end of the simulation, and the average fraction of open RYRs approached the maximum steady-state PO of the RYR model. This was the case in 25% of the 40 tested combinations of CRS determinants—CRS sizes, RYR densities, and iCa values—that produced sparks. About 27% of conditions, characterized by somewhat lower coupling strengths, resulted in partial spark deactivation; however, the largest fraction of CRS combinations (48%) led to sparks fully terminating within 100 ms. It should be noted that, at the lowest coupling strengths, sparks did not occur.

Myocyte contraction depends on the total amount of calcium ions released at the CRSs. Fig. 8 B shows the dependence of the amount of released calcium during CREs on the coupling strength φ. For quarks and blips, this amount, nCa, decreased with increasing φ, with no dependence on the number of RYRs in CRS. In the case of sparks, nCa increased with increasing φ, but even more steeply for larger NRYR. These differences are because, at low coupling strength, the activation of blips shortens the duration of quarks. At large coupling strengths, increased activation of a spark results in shorter durations of blips. It is worth noting that the amount of calcium released in quarks, blips, and sparks differs by orders of magnitude (Fig. 8 B).

The amount of released calcium predicted by the inspected CRS models can be compared with experimental estimates. In the case of a narrow CRS with 20 RYRs and iCa of 0.25 pA, used as an example in the above analyses, the amount of released Ca2+ was 0.002, 0.02, and 0.3 amol for the average quark, blip, and spark, respectively. An average CRE released 0.1 amol Ca2+. These numbers can be compared with the amount of 0.1 amol necessary to increase the free cytosolic calcium from 100 nM to 1 µM in a volume served by one couplon (100 µM total Ca2+ in 1 µm3; Bers, 2001; Soeller et al., 2007), and with the 0.06 amol per spark, estimated from Janíček et al. (2012). For comparison, the amount of calcium stored in one cisterna can be estimated as 0.07 amol based on the published data (Hake et al., 2012).

### The effectiveness of spark initiation depends on the initiating RYR

In simulation studies on contact network CRSs, Walker et al., 2014, 2015 showed that the efficiency of spark activation depends on the position of the initiating RYR. A substantial increase in spark fidelity was observed only if the initiating RYR channel was located in the core of the CRS, since RYRs at the periphery had fewer neighbor RYRs than RYRs in the core of the CRS. Although a similar phenomenon could also be expected in less dense CRSs, the extent of RYR recruitment in CRSs of different internal disposition could not be predicted; therefore, we analyzed the recruitment of RYRs by RYR openings at all RYR positions in the CRS for the defined range of CRS determinants (NRYR, AC/S, and iCa). The initiating RYR, characterized by its coupling strength φi, affected distinctly the relative occurrence and the average relative amplitude of CRS responses $(A¯rel),$ as summarized in Fig. 9. It is apparent that differences between RYRs with higher and lower levels of local clustering (i.e., with higher and lower vicinity) become more apparent in CRSs with lower group RYR vicinity (insets in individual panels of Fig. 9).

A typical compact, narrow, and wide CRS with NRYR of 20, 40, and 80 RYRs per CRS was analyzed for the whole range of studied iCa values (0.04, 0.15, 0.25, and 0.4 pA). Altogether, 560 RYRs in 12 CRSs were examined. The dependences of CRE characteristics on φi are shown in Fig. 10. The average relative amplitude $A¯reli$ and the fraction of sparks $Fsi$ increased, while the fraction of quarks $Fqi$ decreased with increasing φi. The fractional occurrence of blips $Fbi$ showed a maximum of ∼0.4 in the lower intermediate range of φi values. The dependences of CRE characteristics on φi were approximated by the respective functions derived from Eqs. 7 and 9. The value of the exponent α was found by relating all observed event characteristics of all simulated CRSs to the variable φi in parallel. The best-fit parameters are summarized in Table 3. The RYR coupling strength φi allowed us to line up relationships across spark characteristics and CRS density types at different NRYR and iCa.

### CRS spark activity results from synergy between RYR vicinity and iCa

The previous sections demonstrated a functional relationship between spark activity and the coupling strength characterizing individual RYRs (φi) as well as whole CRSs (φ). Since φ values are products of two quantities (Eqs. 7 and 8), we need to know the contribution of each to the formation of CREs. Fig. 10 summarizes the calculated dependences of relative CRE amplitudes and relative occurrences of quarks, blips, and sparks on a wide range of iCa and v values. For comparison, the vicinities of CRS models used for the simulations are marked on the right side. For all CRE characteristics, the hyperbolic relationships draw limits in which the CRSs of a specific density type can be effectively active. These relationships also indicate under what geometric conditions CRSs can be effectively controlled by calcium flux, as well as under what conditions CRSs would not activate at any physiologically possible iCa. The amplification of the effect of v by iCa, and vice versa, is a manifestation of the synergy between iCa and v.

### Calcium fluorescence signals calculated from CREs

The question is, how are the simulated CREs related to calcium sparks observed in cardiac myocytes by confocal fluorescence microscopy. The fluorescence signal of calcium sparks is proportional to local [Ca2+]. Since it is quantitatively related to calcium release flux, it can be calculated from the simulated CREs, e.g., as described by Smith et al. (1998) and Zahradníková et al. (2007). Fluorescence intensities calculated for the averaged $NO(t)$ records of the quark, blip, and spark (red traces in Fig. 8 A, left) are shown in Fig. 11 as x-t images in confocal microscopy experiments. The characteristics of the respective fluorescence signals are summarized in Table 4. The fluorescence of an average quark is below the standard resolution limit of confocal microscopy, as expected (Lipp and Niggli, 1996; Niggli and Egger, 2002). According to our estimates, a quark could be visible if the calcium current were large and the RYR openings were simultaneously long, or if the cytosolic buffering capacity for calcium were small, or if the free diffusion space were restricted (e.g., by membranous organelles). The fluorescence of an average blip resembles quarky calcium releases, as observed by several laboratories (Brochet et al., 2011; Shang et al., 2014). The fluorescence of an average spark resembles, by its appearance, typical calcium sparks, as observed by many laboratories. For comparison, Zahradníková et al. (2007) reported spark amplitude of 1.75 ΔF/F0 units, a TTP of 16 ms, and full duration at half-maximun of 38 ms. Thus, these calculations confirmed the adequacy of the use of simulated CREs for inferences on CRS behavior and mechanisms of local calcium release.

### Local calcium release is shaped by RYR patterns and CRS arrangement

The finding by Asghari et al. (2014) of a direct relationship between RYR arrangement in dyads and the concentration of Mg2+ ions opened the question of its role in the physiological regulation of excitation–contraction coupling in cardiac myocytes. We tested whether the arrangement of RYRs into ChB and SbS patterns (see Table 1) would generate physiologically meaningful differences in CRS activity in the given simulation settings. For narrow CRSs and iCa of 0.25 pA, the relative occurrences of quarks, blips, and sparks showed only minor differences between tested patterns when NRYR was the same (Table 5). Although statistical testing indicated that Fq, Fb, or Fs values were significantly different between some CRS patterns at some NRYR values, these differences were not consistent for CRSs with different NRYR.

Generally, the differences in the fraction of quarks, blips, and sparks between CRSs of different pattern types, but of the same NRYR, were much smaller than those between CRSs of the same pattern type but with different NRYR. These results could be expected since the only interaction factor controlling the activation of RYRs implemented in this study was the local calcium concentration. Despite the steep gradient at an open RYR, the concentration of Ca2+ was similar at the nearest neighbor RYR whether in SbS or the slightly more distant ChB position. Nevertheless, the coupled gating or another type of RYR–RYR interaction might contribute differently to the activation of CRSs with different pattern types. To reveal this, an extended RYR gating model would need to be inspected.

In previous sections, we analyzed the effects of RYR arrangement at the level of a single CRS. In real experiments, the spatial resolution of calcium sparks is less than the physical size of a CRS; thus, the proximate CRSs cannot be resolved. The question is whether assembling RYRs into either one or a few proximate, but independent, CRSs could be distinguished by their activity pattern. We compared the formation of sparks in three scenarios differing in the arrangement of 40 RYRs into one, two, or four independent CRSs (Table 6) so that no CRS is affected by Ca2+ diffusion from the activated CRSs (nota bene: in multicluster models of dyads, the approximation of independent RYR clusters would correspond to edge-to-edge separation of RYR clusters larger than 150 nm [Macquaide et al., 2015; Kolstad et al., 2018]). The single large, wide CRS with 40 RYRs represents the first scenario, in which all RYRs are coupled by diffusing calcium. The second scenario comprises two independent nonidentical CRSs, each of the narrow type and with 20 RYRs. The third scenario comprises four independent nonidentical CRSs, each with 10 RYRs in the compact arrangement. In all scenarios, the simulations were started by activation of an RYR in only one of the CRSs. This CRE activation procedure tallies the frequency of initializing RYR opening per scenario since the total number of RYRs is the same in each scenario.

Changes in CRS activity calculated for these scenarios showed dramatic differences in the averaged CRE responses (Table 6). The single large CRS with 40 RYRs of low vicinity produced sparks rather rarely, just in five of 100 initiating RYR openings. The two 20 RYR CRSs of the narrow type (i.e., with a higher RYR vicinity) produced 29 sparks, and the four 10 RYR CRSs of the compact type—and thus of the highest RYR vicinity—produced as many as 60 sparks of 100 initiating RYR openings. In these scenarios, the average number of simultaneously open RYRs in sparks, $N¯Opeak$ (Table 6), and the amount of calcium released per spark, nCa,s (not shown), were predictably reduced following the reduction of NRYR per one CRS. However, since the decrease in the spark amplitude was counteracted by the increase in the spark frequency, the amount of calcium released by sparks per one initiating RYR opening (Fs × nCa,s; Table 6) was substantially larger in the two-CRS case relative to the single-CRS case, but it was intermediate in the four-CRS case. On the contrary, compacting the RYRs into a few small CRSs resulted in a smaller amount of sparkless release (Fq × nCa,q, and Fb × nCa,b; Table 6).

Analysis of RYR activity in individual CRSs shown in Table 6 revealed that the propensity of RYRs to initiate a spark, $Fsi,$ increased with increasing RYR vicinity, in agreement with previous simulations (Fig. 9). The fluorescence of averaged sparks (Table 6) corresponded to the time course of RYR activity in the respective CRSs (Fig. 8). Compared with the one-CRS scenario, the two- and four-CRS scenarios produced sparks with a progressively smaller amplitude; however, the spark produced by the two-CRS scenario was longer, while that produced by the four-CRS scenario was shorter. These changes are in line with previous findings on the relationships between spark amplitude and duration, on one hand, and the vicinity and the number of interacting RYRs, on the other hand (Fig. 8).

Altogether, these simulations suggest that the frequency and amplitude of calcium sparks generated by a set of unresolved CRSs depend strongly on the number, size, and density type of these CRSs. Moreover, the fine structure of CRSs determines the ratio of “invisible” nonspark calcium release (Ca leak) and “visible” spark-related calcium release.

Recent experimental and theoretical studies indicated the importance of dyad structure and RYR clustering on spontaneous calcium spark formation (Cannell et al., 2013; Walker et al., 2014, 2015; Kolstad et al., 2018; Cosi et al., 2019) as well as on triggered calcium release (Zahradnik et al., 2013; Cosi et al., 2019; Novotová et al., 2020). In this study, we used a large set of CRS geometries that approximated the structural data by constrained random distribution of RYRs in CRS (Iaparov et al., 2017; Iaparov et al., 2019).

This allowed us to solve the role of geometric factors of cardiac dyads in the formation of elementary CREs. The modeling approach allowed for wide-ranging exploration of dyad-defining determinants implemented into numerous models of CRSs, which would not be feasible in real experiments. The scope was limited to a realistic range of determinants, as estimated thus far. The in silico experiments emulated spontaneous CREs observed in resting cardiac myocytes. The openings of a single RYR in the CRS generated CREs that were recorded and analyzed. Under a majority of CRS conditions, three types of CREs, namely, quarks, blips, and sparks, were generated. Their relative occurrence, amplitude, and kinetics revealed high cooperativity between the number, arrangement, and single-channel calcium current of RYRs in the formation of CRS responses. This finding was explained by introducing the RYR coupling strength φ, defined as a weighted contribution of microscopic CRS determinants—the single-channel calcium current and the vicinity of RYRs in the CRS, where the RYR vicinity is a new determinant that characterizes the arrangement of RYRs in the CRS, invariant to CRS size.

RYR vicinity (Eqs. 4 and 5) characterizes distances between RYRs in the group by a single numerical value. Unlike the RYR surface density or the number of RYRs in the CRS, the vicinity resolves the differences in the distribution of RYRs in dyads. In other words, two CRSs that have the same number of RYRs and the same surface density may have different group vicinity values. Likewise, two RYRs in two CRSs of different surface densities may have the same RYR vicinity value. To generalize, the vicinity is a measure of the RYR interaction potential since it increases with increasing RYR clustering. The previously used concept of adjacency or nearest neighborhood (Walker et al., 2014) attributes the interaction potential only to close neighbors, resulting in discretized levels of RYR–RYR interactions. In contrast to the RYR surface density or mean RYR occupancy used in previous models (Cannell et al., 2013; Cosi et al., 2019), the vicinity characterizes the whole CRS (Eq. 6) as well as individual RYRs (Eq. 5) according to their global or local geometric dispositions.

We defined the coupling strength φ (Eqs. 7 and 8) as an internal determinant of CRS productivity, which combines the vicinity of interacting elements (RYRs) with the strength of the interaction agents ([Ca2+] and weighted iCa). It turned out that the amplitudes, kinetic characteristics, and probabilities of occurrence of all CRE types correlated well with this CRS determinant. Importantly, the optimal correlation was obtained with a single value of the weight factor α (Figs. 7 B and 8). Therefore, coupling strength φ could be thought of as a measure of effective calcium concentration in the dyadic gap (Eq. 9), although its dimension Aα m−1 corresponds to a linear current density or calcium flux per distance or perimeter. What is important, however, is that the same relationship could be applied to responses of CRS models built on the experimentally determined (Walker et al., 2015) and algorithm-generated RYR distributions (Iaparov et al., 2019) at a wide range of RYR surface densities. Taken together, these multiparametric correlations on the grounds of generalized RYR geometries support the intuitive understanding that the vicinity and iCa of RYRs determine the activity of CRS in synergy (Fig. 10).

The RYR coupling strength φ has some similarity to the previously introduced concept of coupling strength c* (DeRemigio and Smith, 2005; Groff and Smith, 2008), which was used as an ad hoc input parameter in calculations to characterize the average calcium concentration increase in the CRS caused by the opening of a single channel. Therefore, this parameter cannot be decomposed into a geometric and a functional component. The coupling strength c* underestimates the interaction of RYRs in noncompact CRS (Iaparov et al., 2019), and thus cannot consolidate the correlations between the properties of CREs and the determinants of general CRSs. In contrast, the RYR coupling strength φ conveniently integrates the geometric and functional aspects of CRS activity. Its effective value is determined ex post by approximating the characteristics of simulated CREs by a function of RYR vicinity, calcium current, and the weighing factor α. The coupling strengths φ and φi deliver the synergic character of RYR clustering and calcium current and thus allow for prediction of the relative occurrence of quarks, blips, and sparks, as well as their amplitudes and activation kinetics in the whole range of studied CRS characteristics; however, individual responses of a CRS may vary broadly, as does the chance of RYR opening. Taken from the other side, CRS and RYR coupling strength determine the propensity of CRSs to generate specific types of CREs, their amplitudes, and their time courses. They may predict whether a specific CRS will generate preferentially quarks, blips, or sparks (Fig. 10). The coupling strength makes it possible to explain the deterministic relationship between the stochastic processes of RYR gating, the random RYR distribution, and calcium diffusion, on the one hand, and the generation of CREs, on the other hand. To fully accomplish these goals, we need new models of CRS with an improved definition of the control of RYR gating by luminal calcium at both luminal (Terentyev et al., 2002; Laver, 2007; Tencerová et al., 2012) and cytosolic sites (Sobie et al., 2002; Laver, 2007) that would account for the fact that the coupling strengths φ and φi decline with SR depletion. Another open question is the existence and relevance of direct RYR–RYR interactions that occur during coupled gating (Stern et al., 1999; Marx et al., 2001) or other yet unrevealed coupling mechanisms.

All CRS models explored in this study responded to activation by CREs of amplitudes and time courses that were highly variable. The CRE variability was high in CRSs of any RYR arrangement, despite the single-channel current value being constant during CREs, and thus could be fully attributed to the gating and group behavior of RYRs. However, the geometric disposition of RYRs substantially modulated the characteristics of CREs. CREs of all three types were generated by any CRS model, although not at any iCa value; however, their relative occurrences were widely different. Moreover, repeated activation of the same CRS by the same initiating RYR also produced all three types of CREs (Figs. 4 and 6). Therefore, we suppose that the various types of local CREs observed in cardiac myocytes under otherwise the same conditions represent a fundamental variability of CRS responses. In other words, the experiments indicating a broad range of calcium spark amplitudes (Wang et al., 2004; Zima et al., 2008; Brochet et al., 2011; Janíček et al., 2012) and durations (Zima et al., 2008; Brochet et al., 2011), or even opening of one or a few RYRs (Wang et al., 2004; Janíček et al., 2012), have support in the theory of CRS operation.

The theoretical study by Stern et al. (1999) predicted that calcium release may terminate even in the absence of a specific termination mechanism due to stochastic attrition. Recent studies of CRS at a constant RYR current predicted a threshold value of calcium current below which any CRS event will terminate (Maltsev et al., 2019; Gillespie, 2020). This study revealed the stochastic termination dominating in blips and prominent in sparks. The picture of calcium spark termination, emerging from our simulations, points to a complex interplay between CRS determinants and kinetics of RYR gating. Importantly, however, sparks produced by the tested CRSs differed too broadly in their duration. This is at odds with the observed characteristics of spark termination. Without luminal regulation, typical calcium sparks in large dyads and at high SR calcium loads would be much higher and would last much longer than physiologically useful. In myocytes, therefore, spark termination has to involve control mechanisms that would set properly the duration of local CREs.

### Experimental implications

The results of these simulations may be of help to settle some experimental problems. For instance, they indicate that the single openings of RYRs, classified here as quarks, produce fluorescence signals too small to be detectable as local CREs in myocytes under conditions close to physiological. The activity of a small group of RYRs produces blips that release calcium ions in amounts sufficient to generate calcium signals that are close to the detection limit. The spark-like activity of a larger group of RYRs yields large calcium signals under a broad range of CRS settings. In noncompact CRSs, the rate of RYR recruitment seems to be relatively slow, and thus well-designed fluorescence measurements may characterize the kinetics of calcium sparks even when convoluted with the kinetics of calcium binding by a fluorescent dye. This may help to better resolve the fine structure of the CRE source, since the amplitude and kinetics of calcium sparks strongly depend on the distribution of RYRs within one or several CRSs (Table 6).

According to these simulations, brief sparks are produced by CRSs with either a small number of RYRs or a low group vicinity or a smaller RYR calcium current. Slowly terminating sparks are produced by CRSs with either a large number of RYRs or a high group vicinity or a high iCa. Large sparks are produced by CRSs with a large number of RYRs, and their amplitude can be further increased by a large coupling strength (Fig. 8).

It is well documented that, during a calcium spark, a substantial local depletion of luminal Ca2+ concentration occurs (Brochet et al., 2005), amounting to 40–90% (Zima et al., 2008; Kong et al., 2013). Considering the example of narrow CRS with 20 RYRs and iCa of 0.25 pA, an ∼40% decrease in iCa would result in transformation of 72% of sparks to blips or quarks. According to Kong et al. (2013), 40% depletion may be achieved 10 ms after the start of a CRE. In the example CRS, however, 90% of peak activation is achieved in <5 ms, and therefore even such a rapid depletion would affect only slightly the number of recruited RYRs or the fraction of sparks, but would significantly curtail the duration of calcium release flux. Similar considerations propose that dyadic SR cisternae with low Ca content or with lower RYR Ca conductance may have a much lower probability of calcium spark activation and, conversely, a higher contribution of nonspark calcium release in the form of quarks and blips, as previously observed by Zima et al. (2010).

For calcium release experiments, especially those aimed at pathological mechanisms, of consideration can be the changes in RYR distribution. According to these simulations, the change in RYR vicinity affects the CRS propensity to generate sparks (Table 6), which can be easily measured and could be correctly interpreted if the SR calcium content was constant. The RYR vicinity also correlates with the rate of spark activation (Fig. 8 A) since the coupling strength is directly proportional to the vicinity. This is in agreement with previous experimental data and simulations (Kolstad et al., 2018).

The calculated amount of calcium released per quark and per blip is 10–100 times smaller than that released per spark according to these simulations; however, it may contribute significantly to calcium leak from the SR when their occurrence is high. Importantly, the ratio of invisible nonspark calcium release (Ca leak) and visible spark-related calcium release depends strongly on the distribution of a given number of RYRs into one or several optically unresolvable CRSs (Table 6).

The results presented in Fig. 8 may partially explain the effect of the store overload-induced calcium release (Kong et al., 2007), since they predict a strong increase in spark occurrence and duration upon an increase of iCa. Similar considerations may be applied to observations and interpretation of calcium waves in cardiac myocytes that accompany the development of arrhythmias in injured or failing myocardia.

David A. Eisner served as editor.

The study was performed using the Uran supercomputer of the Krasovskii Institute of Mathematics and Mechanics.

The research was supported by Russian Foundation for Basic Research according to research project no. 18-31-00153 (development of the computational model of CRS); by the Government of the Russian Federation (program 02.A03.21.0006); by the Ministry of Science and Higher Education of the Russian Federation (project no. FEUZ-2020-0054); by the Slovak Research and Development Agency (project APVV-15-0302); by SAV-TUBITAK (project JRP/2019/836/RyRinHeart); and by the Operational Program Integrated Infrastructure for the project “Long-term strategic research of prevention, intervention, and mechanisms of obesity and its comorbidities,” ITMS: 313011V344, co-financed by the European Regional Development Fund.

The authors declare no competing financial interests.

Author contributions: B.I. Iaparov: conceptualization, methodology, software, data curation, validation, formal analysis, investigation, writing – original draft, writing – review and editing, project administration, and funding acquisition. I. Zahradník: conceptualization, methodology, validation, investigation, writing – original draft, writing – review and editing, and visualization. A.S. Moskvin: conceptualization, supervision, and writing – review and editing. A. Zahradníková: conceptualization, methodology, software, data curation, validation, formal analysis, investigation, writing – original draft, writing – review and editing, visualization, supervision, project administration, and funding acquisition.

Asghari
,
P.
,
D.R.
Scriven
,
S.
Sanatani
,
S.K.
Gandhi
,
A.I.
Campbell
, and
E.D.
Moore
.
2014
.
Nonuniform and variable arrangements of ryanodine receptors within mammalian ventricular couplons
.
Circ. Res.
115
:
252
262
.
Asghari
,
P.
,
D.R.
Scriven
,
M.
Ng
,
P.
Panwar
,
K.C.
Chou
,
F.
van Petegem
, and
E.D.
Moore
.
2020
.
Cardiac ryanodine receptor distribution is dynamic and changed by auxiliary proteins and post-translational modification
.
eLife.
9
:e51602.
,
D.
,
I.D.
Jayasinghe
,
L.
Lam
,
S.
Rossberger
,
M.B.
Cannell
, and
C.
Soeller
.
2009
.
Optical single-channel resolution imaging of the ryanodine receptor distribution in rat cardiac myocytes
.
106
:
22275
22280
.
Bers
,
D.M.
2001
.
Sources and sinks of activator calcium
. In
Excitation-contraction coupling and cardiac contractile force.
,
Dordrecht
.
39
62
.
Brochet
,
D.X.
,
D.
Yang
,
A.
Di Maio
,
W.J.
Lederer
,
C.
Franzini-Armstrong
, and
H.
Cheng
.
2005
.
Ca2+ blinks: rapid nanoscopic store calcium signaling
.
102
:
3099
3104
.
Brochet
,
D.X.
,
W.
Xie
,
D.
Yang
,
H.
Cheng
, and
W.J.
Lederer
.
2011
.
Quarky calcium release in the heart
.
Circ. Res.
108
:
210
218
.
Cabra
,
V.
,
T.
Murayama
, and
M.
Samsó
.
2016
.
Ultrastructural analysis of self-associated RyR2s
.
Biophys. J
.
110
:
2651
2662
.
Cannell
,
M.B.
,
H.
Cheng
, and
W.J.
Lederer
.
1995
.
The control of calcium release in heart muscle
.
Science.
268
:
1045
1049
.
Cannell
,
M.B.
,
C.H.
Kong
,
M.S.
Imtiaz
, and
D.R.
Laver
.
2013
.
Control of sarcoplasmic reticulum Ca2+ release by stochastic RyR gating within a 3D model of the cardiac dyad and importance of induction decay for CICR termination
.
Biophys. J.
104
:
2149
2159
.
Cheng
,
H.
,
W.J.
Lederer
, and
M.B.
Cannell
.
1993
.
Calcium sparks: elementary events underlying excitation-contraction coupling in heart muscle
.
Science.
262
:
740
744
.
Cheng
,
H.
,
M.R.
Lederer
,
W.J.
Lederer
, and
M.B.
Cannell
.
1996
.
Calcium sparks and [Ca2+]i waves in cardiac myocytes
.
Am. J. Physiol.
270
:
C148
C159
.
Cleemann
,
L.
,
W.
Wang
, and
M.
.
1998
.
Two-dimensional confocal images of organization, density, and gating of focal Ca2+ release sites in rat cardiac myocytes
.
95
:
10984
10989
.
Cosi
,
F.G.
,
W.
Giese
,
W.
Neubert
,
S.
Luther
,
N.
Chamakuri
,
U.
Parlitz
, and
M.
Falcke
.
2019
.
Multiscale modeling of dyadic structure-function relation in ventricular cardiac myocytes
.
Biophys. J
.
117
:
2409
2419
.
Crank
,
J.
1975
.
Mathematics of Diffusion.
Oxford University Press
,
Oxford, UK
. 414 pp.
DeRemigio
,
H.
, and
G.D.
Smith
.
2005
.
The dynamics of stochastic attrition viewed as an absorption time on a terminating Markov chain
.
Cell Calcium.
38
:
73
86
.
Drum
,
B.M.
,
C.
Yuan
,
A.
de la Mata
,
N.
Grainger
, and
L.F.
Santana
.
2020
.
Junctional sarcoplasmic reticulum motility in adult mouse ventricular myocytes
.
Am. J. Physiol. Cell Physiol.
318
:
C598
C604
.
Eisner
,
D.A.
,
J.L.
Caldwell
,
K.
Kistamás
, and
A.W.
Trafford
.
2017
.
Calcium and excitation-contraction coupling in the heart
.
Circ. Res
.
121
:
181
195
.
Fabiato
,
A.
1983
.
Calcium-induced release of calcium from the cardiac sarcoplasmic reticulum
.
Am. J. Physiol.
245
:
C1
C14
.
Franzini-Armstrong
,
C.
,
F.
Protasi
, and
V.
Ramesh
.
1998
.
Comparative ultrastructure of Ca2+ release units in skeletal and cardiac muscle
.
853
(
1 CARDIAC SARCO
):
20
30
.
Franzini-Armstrong
,
C.
,
F.
Protasi
, and
V.
Ramesh
.
1999
.
Shape, size, and distribution of Ca2+ release units and couplons in skeletal and cardiac muscles
.
Biophys. J
.
77
:
1528
1539
.
Galice
,
S.
,
Y.
Xie
,
Y.
Yang
,
D.
Sato
, and
D.M.
Bers
.
2018
.
Size matters: Ryanodine receptor cluster size affects arrhythmogenic sarcoplasmic reticulum calcium release
.
J. Am. Heart Assoc
.
7
. e008724.
Gillespie
,
D.
2020
.
Recruiting RyRs to open in a Ca2+ release unit: Single-RyR gating properties make RyR group dynamics
.
Biophys J
.
118
:
232
242
.
Gillespie
,
D.T.
1977
.
Exact stochastic simulation of coupled chemical reactions
.
J. Phys. Chem
.
81
:
2340
2361
.
Groff
,
J.R.
, and
G.D.
Smith
.
2008
.
Ryanodine receptor allosteric coupling and the dynamics of calcium sparks
.
Biophys. J.
95
:
135
154
.
Györke
,
I.
, and
S.
Györke
.
1998
.
Regulation of the cardiac ryanodine receptor channel by luminal Ca2+ involves luminal Ca2+ sensing sites
.
Biophys. J.
75
:
2801
2810
.
Hake
,
J.
,
A.G.
Edwards
,
Z.
Yu
,
P.M.
Kekenes-Huskey
,
A.P.
Michailova
,
J.A.
McCammon
,
M.J.
Holst
,
M.
Hoshijima
, and
A.D.
McCulloch
.
2012
.
Modelling cardiac calcium sparks in a three-dimensional reconstruction of a calcium release unit
.
J. Physiol.
590
:
4403
4422
.
Hayashi
,
T.
,
M.E.
Martone
,
Z.
Yu
,
A.
Thor
,
M.
Doi
,
M.J.
Holst
,
M.H.
Ellisman
, and
M.
Hoshijima
.
2009
.
Three-dimensional electron microscopy reveals new details of membrane systems for Ca2+ signaling in the heart
.
J. Cell Sci.
122
:
1005
1013
.
Hiess
,
F.
,
A.
Vallmitjana
,
R.
Wang
,
H.
Cheng
,
H.E.
ter Keurs
,
J.
Chen
,
L.
,
R.
Benitez
, and
S.R.
Chen
.
2015
.
Distribution and function of cardiac ryanodine receptor clusters in live ventricular myocytes
.
J. Biol. Chem
.
290
:
20477
20487
.
Hou
,
Y.
,
I.
Jayasinghe
,
D.J.
Crossman
,
D.
, and
C.
Soeller
.
2015
.
Nanoscale analysis of ryanodine receptor clusters in dyadic couplings of rat cardiac myocytes
.
J. Mol. Cell. Cardiol.
80
:
45
55
.
Iaparov
,
B.I.
,
S.Y.
Khamzin
,
A.S.
Moskvin
, and
O.E.
Solovyova
.
2017
.
Mathematical modeling shows the frequency of Ca2+ sparks in cells depends on the ryanodine receptor’s arrangement
.
Procedia Comput. Sci
.
119
:
190
196
.
Iaparov
,
B.I.
,
A.S.
Moskvin
,
I.
, and
A.
.
2019
.
Stochastic and deterministic approaches to modelling calcium release in cardiac myocytes at different spatial arrangements of ryanodine receptors
.
Eur. Biophys. J.
48
:
579
584
.
Iaparov
,
B.I.
,
I.
,
A.S.
Moskvin
, and
A.
.
2020
.
Synergy of calcium release site determinants in control of calcium release events in cardiac myocytes
.
bioRxiv
. doi: . (Preprint posted August 27, 2020).
Janíček
,
R.
,
A.
Jr
.,
E.
Poláková
,
J.
Pavelková
,
I.
, and
A.
.
2012
.
Calcium spike variability in cardiac myocytes results from activation of small cohorts of ryanodine receptor 2 channels
.
J. Physiol.
590
:
5091
5106
.
Jayasinghe
,
I.
,
A.H.
Clowsley
,
R.
Lin
,
T.
Lutz
,
C.
Harrison
,
E.
Green
,
D.
,
L.
Di Michele
, and
C.
Soeller
.
2018
.
True molecular scale visualization of variable clustering properties of ryanodine receptors
.
Cell Rep
.
22
:
557
567
.
,
T.R.
,
J.
van den Brink
,
N.
MacQuaide
,
P.K.
Lunde
,
M.
Frisk
,
J.M.
Aronsen
,
E.S.
Norden
,
A.
Cataliotti
,
I.
,
O.M.
Sejersted
, et al
.
2018
.
Ryanodine receptor dispersion disrupts Ca2+ release in failing cardiac myocytes
.
eLife.
7
:e39427.
Kong
,
H.
,
R.
Wang
,
W.
Chen
,
L.
Zhang
,
K.
Chen
,
Y.
Shimoni
,
H.J.
Duff
, and
S.R.
Chen
.
2007
.
Skeletal and cardiac ryanodine receptors exhibit different responses to Ca2+ overload and luminal ca2+
.
Biophys. J.
92
:
2757
2770
.
Kong
,
C.H.
,
D.R.
Laver
, and
M.B.
Cannell
.
2013
.
Extraction of sub-microscopic Ca fluxes from blurred and noisy fluorescent indicator images with a detailed model fitting approach
.
PLOS Comput. Biol.
9
:e1002931.
Laver
,
D.R.
2007
.
Ca2+ stores regulate ryanodine receptor Ca2+ release channels via luminal and cytosolic Ca2+ sites
.
Biophys. J.
92
:
3541
3555
.
Laver
,
D.R.
, and
B.N.
Honen
.
2008
.
Luminal Mg2+, a key factor controlling RYR2-mediated Ca2+ release: cytoplasmic and luminal regulation modeled in a tetrameric channel
.
J. Gen. Physiol.
132
:
429
446
.
Laver
,
D.R.
,
T.M.
Baynes
, and
A.F.
Dulhunty
.
1997
.
Magnesium inhibition of ryanodine-receptor calcium channels: evidence for two independent mechanisms
.
J. Membr. Biol.
156
:
213
229
.
Lipp
,
P.
, and
E.
Niggli
.
1996
.
Submicroscopic calcium signals as fundamental events of excitation--contraction coupling in guinea-pig cardiac myocytes
.
J. Physiol.
492
:
31
38
.
López-López
,
J.R.
,
P.S.
Shacklock
,
C.W.
Balke
, and
W.G.
Wier
.
1994
.
Local, stochastic release of Ca2+ in voltage-clamped rat heart cells: visualization with confocal microscopy
.
J. Physiol.
480
:
21
29
.
López-López
,
J.R.
,
P.S.
Shacklock
,
C.W.
Balke
, and
W.G.
Wier
.
1995
.
Local calcium transients triggered by single L-type calcium channel currents in cardiac cells
.
Science.
268
:
1042
1045
.
Louch
,
W.E.
,
J.T.
Koivumäki
, and
P.
Tavi
.
2015
.
Calcium signalling in developing cardiomyocytes: implications for model systems and disease
.
J. Physiol.
593
:
1047
1063
.
Lukyanenko
,
V.
, and
S.
Gyorke
.
1999
.
Ca2+ sparks and Ca2+ waves in saponin-permeabilized rat ventricular myocytes
.
J. Physiol.
521
:
575
585
.
Macková
,
K.
,
A.
Jr
.,
M.
Hoťka
,
B.
Hoffmannová
,
I.
, and
A.
.
2017
.
Calcium release-dependent inactivation precedes formation of the tubular system in developing rat cardiac myocytes
.
Eur. Biophys. J.
46
:
691
703
.
Macquaide
,
N.
,
H.T.
Tuan
,
J.
Hotta
,
W.
Sempels
,
I.
Lenaerts
,
P.
Holemans
,
J.
Hofkens
,
M.S.
Jafri
,
R.
Willems
, and
K.R.
Sipido
.
2015
.
Ryanodine receptor cluster fragmentation and redistribution in persistent atrial fibrillation enhance calcium release
.
Cardiovasc. Res.
108
:
387
398
.
Maltsev
,
A.V.
,
M.D.
Stern
, and
V.A.
Maltsev
.
2019
.
Mechanisms of calcium leak from cardiac sarcoplasmic reticulum revealed by statistical mechanics
.
Biophys J
.
116
:
P2212
P2223
.
Marx
,
S.O.
,
J.
Gaburjakova
,
M.
Gaburjakova
,
C.
Henrikson
,
K.
Ondrias
, and
A.R.
Marks
.
2001
.
Coupled gating between cardiac calcium release channels (ryanodine receptors)
.
Circ. Res.
88
:
1151
1158
.
Matveev
,
V.
,
A.
Sherman
, and
R.S.
Zucker
.
2002
.
New and corrected simulations of synaptic facilitation
.
Biophys. J.
83
:
1368
1373
.
Meissner
,
G.
2002
.
Regulation of mammalian ryanodine receptors
.
Front. Biosci.
7
:
d2072
d2080
.
Meissner
,
G.
2004
.
Molecular regulation of cardiac ryanodine receptor ion channel
.
Cell Calcium.
35
:
621
628
.
Naraghi
,
M.
, and
E.
Neher
.
1997
.
Linearized buffered Ca2+ diffusion in microdomains and its implications for calculation of [Ca2+] at the mouth of a calcium channel
.
J. Neurosci.
17
:
6961
6973
.
Niggli
,
E.
, and
M.
Egger
.
2002
.
Calcium quarks
.
Front. Biosci.
7
:
d1288
d1297
.
Novotova
,
M.
,
A.
, and
I.
.
2009
.
Remodelling of E-C coupling units induced by a single dose of isoproterenol in rats
.
J. Physiol. Sci.
59
:
313
.
Novotova
,
M.
,
Z.
Nichtova
, and
I.
.
2013
.
Changes in dyad morphogenesis in cardiac myocytes due to isoproterenol-induced myocardial injury
.
Eur. J. Heart Fail.
12
:
S43
. Abstr P1744.
Novotová
,
M.
,
A.
Jr
.,
Z.
Nichtová
,
R.
Kováč
,
E.
Kráľová
,
T.
Stankovičová
,
A.
, and
I.
.
2020
.
Structural variability of dyads relates to calcium release in rat ventricular myocytes
.
Sci. Rep.
10
:
8076
.
Parker
,
I.
,
W.J.
Zang
, and
W.G.
Wier
.
1996
.
Ca2+ sparks involving multiple Ca2+ release sites along Z-lines in rat heart cells
.
J. Physiol.
497
:
31
38
.
Scriven
,
D.R.
,
P.
Asghari
,
M.N.
Schulson
, and
E.D.
Moore
.
2010
.
Analysis of Cav1.2 and ryanodine receptor clusters in rat ventricular myocytes
.
Biophys. J.
99
:
3923
3929
.
Shang
,
W.
,
F.
Lu
,
T.
Sun
,
J.
Xu
,
L.L.
Li
,
Y.
Wang
,
G.
Wang
,
L.
Chen
,
X.
Wang
,
M.B.
Cannell
, et al
.
2014
.
Imaging Ca2+ nanosparks in heart with a new targeted biosensor
.
Circ. Res.
114
:
412
420
.
Smith
,
G.D.
,
J.E.
Keizer
,
M.D.
Stern
,
W.J.
Lederer
, and
H.
Cheng
.
1998
.
A simple numerical model of calcium spark formation and detection in cardiac myocytes
.
Biophys. J.
75
:
15
32
.
Sobie
,
E.A.
,
K.W.
Dilly
,
J.
dos Santos Cruz
,
W.J.
Lederer
, and
M.S.
Jafri
.
2002
.
Termination of cardiac Ca(2+) sparks: an investigative mathematical model of calcium-induced calcium release
.
Biophys. J.
83
:
59
78
.
Sobie
,
E.A.
,
L.S.
Song
, and
W.J.
Lederer
.
2006
.
Restitution of Ca(2+) release and vulnerability to arrhythmias
.
J. Cardiovasc. Electrophysiol.
17
(
s1
,
Suppl 1
):
S64
S70
.
Soeller
,
C.
, and
M.B.
Cannell
.
1997
.
Numerical simulation of local calcium movements during L-type calcium channel gating in the cardiac diad
.
Biophys. J.
73
:
97
111
.
Soeller
,
C.
, and
M.B.
Cannell
.
2004
.
Analysing cardiac excitation-contraction coupling with mathematical models of local control
.
Prog. Biophys. Mol. Biol.
85
:
141
162
.
Soeller
,
C.
,
D.
Crossman
,
R.
Gilbert
, and
M.B.
Cannell
.
2007
.
Analysis of ryanodine receptor clusters in rat and human cardiac myocytes
.
104
:
14958
14963
.
Soeller
,
C.
,
I.D.
Jayasinghe
,
P.
Li
,
A.V.
Holden
, and
M.B.
Cannell
.
2009
.
Three-dimensional high-resolution imaging of cardiac proteins to construct models of intracellular Ca2+ signalling in rat ventricular myocytes
.
Exp. Physiol.
94
:
496
508
.
Song
,
L.S.
,
J.S.
Sham
,
M.D.
Stern
,
E.G.
Lakatta
, and
H.
Cheng
.
1998
.
Direct measurement of SR release flux by tracking ‘Ca2+ spikes’ in rat cardiac myocytes
.
J. Physiol.
512
:
677
691
.
Song
,
L.S.
,
E.A.
Sobie
,
S.
McCulle
,
W.J.
Lederer
,
C.W.
Balke
, and
H.
Cheng
.
2006
.
Orphaned ryanodine receptors in the failing heart
.
103
:
4305
4310
.
Stern
,
M.D.
,
L.S.
Song
,
H.
Cheng
,
J.S.
Sham
,
H.T.
Yang
,
K.R.
Boheler
, and
E.
Ríos
.
1999
.
Local control models of cardiac excitation-contraction coupling. A possible role for allosteric interactions between ryanodine receptors
.
J. Gen. Physiol.
113
:
469
489
.
Stern
,
M.D.
,
E.
Ríos
, and
V.A.
Maltsev
.
2013
.
Life and death of a cardiac calcium spark
.
J. Gen. Physiol.
142
:
257
274
.
Sun
,
X.H.
,
F.
Protasi
,
M.
Takahashi
,
H.
Takeshima
,
D.G.
Ferguson
, and
C.
Franzini-Armstrong
.
1995
.
Molecular architecture of membranes involved in excitation-contraction coupling of cardiac muscle
.
J. Cell Biol.
129
:
659
671
.
Swillens
,
S.
,
P.
Champeil
,
L.
Combettes
, and
G.
Dupont
.
1998
.
Stochastic simulation of a single inositol 1,4,5-trisphosphate-sensitive Ca2+ channel reveals repetitive openings during ‘blip-like’ Ca2+ transients
.
Cell Calcium.
23
:
291
302
.
Tencerová
,
B.
,
A.
,
J.
Gaburjáková
, and
M.
Gaburjáková
.
2012
.
Luminal Ca2+ controls activation of the cardiac ryanodine receptor by ATP
.
J. Gen. Physiol.
140
:
93
108
.
Terentyev
,
D.
,
S.
Viatchenko-Karpinski
,
H.H.
Valdivia
,
A.L.
Escobar
, and
S.
Györke
.
2002
.
Luminal Ca2+ controls termination and refractory behavior of Ca2+-induced Ca2+ release in cardiac myocytes
.
Circ. Res.
91
:
414
420
.
Valent
,
I.
,
A.
,
J.
Pavelková
, and
I.
.
2007
.
Spatial and temporal Ca2+, Mg2+, and ATP2 dynamics in cardiac dyads during calcium release
.
Biochim. Biophys. Acta.
1768
:
155
166
.
Walker
,
M.A.
,
G.S.B.
Williams
,
T.
Kohl
,
S.E.
Lehnart
,
M.S.
Jafri
,
J.L.
Greenstein
,
W.J.
Lederer
, and
R.L.
Winslow
.
2014
.
Superresolution modeling of calcium release in the heart
.
Biophys. J.
107
:
3018
3029
.
Walker
,
M.A.
,
T.
Kohl
,
S.E.
Lehnart
,
J.L.
Greenstein
,
W.J.
Lederer
, and
R.L.
Winslow
.
2015
.
On the adjacency matrix of RyR2 cluster structures
.
PLOS Comput. Biol
.
11
. e1004521.
Wang
,
S.Q.
,
M.D.
Stern
,
E.
Ríos
, and
H.
Cheng
.
2004
.
The quantal nature of Ca2+ sparks and in situ operation of the ryanodine receptor array in cardiac cells
.
101
:
3979
3984
.
Wescott
,
A.P.
,
M.S.
Jafri
,
W.J.
Lederer
, and
G.S.
Williams
.
2016
.
Ryanodine receptor sensitivity governs the stability and synchrony of local calcium release during cardiac excitation-contraction coupling
.
J. Mol. Cell. Cardiol.
92
:
82
92
.
Williams
,
G.S.
,
M.A.
Huertas
,
E.A.
Sobie
,
M.S.
Jafri
, and
G.D.
Smith
.
2007
.
A probability density approach to modeling local control of calcium-induced calcium release in cardiac myocytes
.
Biophys. J.
92
:
2311
2328
.
Williams
,
G.S.
,
A.C.
Chikando
,
H.T.
Tuan
,
E.A.
Sobie
,
W.J.
Lederer
, and
M.S.
Jafri
.
2011
.
Dynamics of calcium sparks and calcium leak in the heart
.
Biophys. J.
101
:
1287
1296
.
Wu
,
H.D.
,
M.
Xu
,
R.C.
Li
,
L.
Guo
,
Y.S.
Lai
,
S.M.
Xu
,
S.F.
Li
,
Q.L.
,
L.L.
Li
,
H.B.
Zhang
, et al
.
2012
.
Ultrastructural remodelling of Ca(2+) signalling apparatus in failing heart cells
.
Cardiovasc. Res.
95
:
430
438
.
Xu
,
Q.L.
,
G.W.
Zhang
,
C.
Zhao
, and
A.M.
An
.
2011
.
A Robust Adaptive Hybrid Genetic Simulated Annealing Algorithm for the global optimization of multimodal functions
.
IEEE.
.
Yin
,
C.C.
,
L.M.
Blayney
, and
F.A.
Lai
.
2005
.
Physical coupling between ryanodine receptor-calcium release channels
.
J. Mol. Biol.
349
:
538
546
.
,
I.
,
S.
Györke
, and
A.
.
2005
.
Calcium activation of ryanodine receptor channels--reconciling RyR gating models with tetrameric channel structure
.
J. Gen. Physiol.
126
:
515
527
.
,
I.
,
A.
,
M.
Novotova
, and
A.
.
2013
.
Changes of calcium release and of dyadic structure in cardiac myocytes due to isoproterenol-induced myocardial injury
.
Eur. J. Heart Fail.
12
:
S43
. Abstr P1749.
,
A.
,
M.
Dura
,
I.
Györke
,
A.L.
Escobar
,
I.
, and
S.
Györke
.
2003
.
Regulation of dynamic behavior of cardiac ryanodine receptor by Mg2+ under simulated physiological conditions
.
Am. J. Physiol. Cell Physiol.
285
:
C1059
C1070
.
,
A.
Jr
.,
E.
Poláková
,
I.
, and
A.
.
2007
.
Kinetics of calcium spikes in rat cardiac myocytes
.
J. Physiol.
578
:
677
691
.
,
A.
,
M.
Gaburjáková
,
J.H.
Bridge
, and
I.
.
2010
a
.
Challenging quantal calcium signaling in cardiac myocytes
.
J. Gen. Physiol.
136
:
581
583
.
,
A.
,
I.
Valent
, and
I.
.
2010
b
.
Frequency and release flux of calcium sparks in rat cardiac myocytes: a relation to RYR gating
.
J. Gen. Physiol.
136
:
101
116
.
Zima
,
A.V.
,
E.
Picht
,
D.M.
Bers
, and
L.A.
Blatter
.
2008
.
Termination of cardiac Ca2+ sparks: role of intra-SR [Ca2+], release flux, and intra-SR Ca2+ diffusion
.
Circ. Res.
103
:
e105
e115
.
Zima
,
A.V.
,
E.
Bovo
,
D.M.
Bers
, and
L.A.
Blatter
.
2010
.
Ca2+ spark-dependent and -independent sarcoplasmic reticulum Ca2+ leak in normal and failing rabbit ventricular myocytes
.
J. Physiol.
588
:
4743
4757
.
Ziman
,
A.P.
,
N.L.
Gómez-Viquez
,
R.J.
Bloch
, and
W.J.
Lederer
.
2010
.
Excitation-contraction coupling changes during postnatal cardiac development
.
J. Mol. Cell. Cardiol.
48
:
379
386
.

## Author notes

A preliminary version of this paper was posted as a preprint in bioRxiv on August 27, 2020.