Construction of Calcium Release Sites in Cardiac Myocytes

Local character of calcium release in cardiac myocytes, as defined by confocal recordings of calcium sparks, implies independent activation of individual calcium release sites based on ryanodine receptor (RyR) channel recruitment. We constructed virtual calcium release sites (vCRSs) composed of a variable number of RyR channels distributed in clusters in accordance with the experimentally observed cluster size distribution. The vCRSs consisted either of a single virtual calcium release unit (vCRU), in which all clusters shared a common dyadic space, or of multiple virtual calcium release units (CRUs) containing one cluster each and having separate dyadic spaces. We explored the stochastic behavior of vCRSs to understand the activation and recruitment of RyRs during calcium sparks. RyRs were represented by the published allosteric gating model that included regulation by cytosolic Ca2+ and Mg2+. The interaction of Mg2+ with the RyR Ca2+-binding sites and the refractory period of vCRSs were optimized to accord with the experimentally observed calcium dependence of calcium spark frequency. The Mg2+-binding parameters of RyRs that provided the best description of spark frequency depended on the number of RyRs assembled in the vCRSs. Adequate inhibitory effect of Mg2+ on the calcium dependence of RyR open probability was achieved if the vCRSs contained at least three clusters. For the distribution of the number of open RyRs in evoked calcium sparks to correspond to the experimentally observed distribution of spark calcium release fluxes, at least three clusters had to share a common virtual CRU, in which ∼3 RyRs open to form an average spark. These results reconcile the small cluster size and stochastic placement of RyRs in the release sites with the estimates of the amount of RyR protein, volume density of calcium release sites, and the size of calcium release sites in rat cardiac myocytes.


INTRODUCTION
Calcium sparks (Cheng et al., 1993), the elementary calcium signals of excitation-contraction coupling in cardiac myocytes, result from activation of ryanodine receptors (RyRs). RyRs channels are clustered at the junctional membrane of the dyadic sarcoplasmic reticulum, thus forming the calcium releasing units (Sun et al., 1995). Stimulation of cardiac myocytes results in synchronous activation of calcium sparks (Cannell et al., 1994) by a calcium influx through the voltage-dependent calcium channels (DHPRs). Sparks occur also in the absence of stimulation under diastolic conditions (Cheng et al., 1993). Such resting calcium sparks are thought to constitute a significant component of calcium leak from the sarcoplasmic reticulum that regulates SR calcium content (Lukyanenko et al., 2001;Maier et al., 2003).
The relationships between RyR activity and calcium spark properties were partially cleared only recently. Specifically, the average number of RyRs that open during the spark was estimated to 2.7-20 (Bridge et al., 1999;Wang et al., 2004;Shkryl et al., 2012), and the calcium dependence of spark frequency was shown to be in line with the calcium dependence of RyR open probability (Lukyanenko et al., 2000;Gusev and Niggli, 2008;Zahradnikova et al., 2010b). The generally accepted image of calcium release units (CRUs) involved large clusters of tightly packed RyRs (Franzini-Armstrong et al., 1999;Soeller et al., 2007). Recent data indicate that both the RyR packing density and the size of CRUs is much smaller (∼14 RyRs per cluster; Baddeley et al., 2009) and highly variable, with RyRs arranged into clusters and superclusters (Baddeley et al., 2009;Hayashi et al., 2009). Quantitative composition of the superclusters has not been determined. It is not known as well, whether adjacent clusters share the same dyadic space and the same terminal cisterna. The observed amount of RyR protein (Bers and Stiffel, 1993) and the density of calcium release sites (Chen-Izu et al., 2006) indicate about 70 RyRs per release site (Soeller et al., 2007;Asghari et al., 2012).
The properties of calcium sparks should be governed by the single channel properties of RyRs, the size of RyR clusters, and the arrangement of clusters into release sites. Previously we have shown (Zahradnikova et al., 2010b) that the relationship between cytosolic calcium concentration and the frequency of spontaneous calcium sparks can be fully explained by gating of single, independent RyRs. Quantitative modeling of these processes with a uniform cluster size of 182 RyRs provided parameters of RyR gating in a very good accordance with both, the experimentally determined gating properties of the RyRs, and the calcium dependence of the The gating scheme of the RyR channel. According to this allosteric homotetrameric RyR model (Zahradnikova et al., 2010b), the single channel properties of RyR channels are described by Eq. 3 and the parameter set given in Table 2. The light blue squares represent the closed states and the orange circles represent the open states of the channel. The first and the second numeral denote the number of calcium and magnesium ions, respectively, bound to the RyR Ca 2+ activation sites. The black (horizontal) and the blue (vertical) lines indicate the calcium-dependent and the magnesium-dependent channel transitions, respectively. The transitions independent of ion binding/unbinding are depicted as dashed lines. All transitions between channel states are reversible and obey the principle of detailed balance. L stays for the low activity L-mode states, and I stays for the inactivated state. Inactivation of the channel by Mg 2+ ions binding to the low-affinity inhibition site is depicted at the right. frequency of spontaneous calcium sparks. It turned out, however, that changes in the size of the CRUs had a strong effect on the simulated calcium spark frequency (Zahradnikova et al., 2010b). Here we explored how the organization of RyRs into small and highly variable clusters can be brought into accordance with RyR gating.
Calcium release flux during calcium sparks has been shown previously to have quantal amplitudes due to opening of only a few RyRs on electrical stimulation (Wang et al., 2004). Quantal distribution of the spark release flux amplitudes is compatible with the model of CRU based on the RyR activation controlled by Mg 2+ /Ca 2+ competition at the RyR activation sites (Zahradnikova et al., 2010b). According to this CRU model, designed with 182 RyRs per CRU, the Mg 2+ ions occupy 92% of the activation sites on RyR channels and prevent explosive RyR activation by calcium. The low rate of Mg 2+ dissociation primes only a few RyRs in the CRU for opening during the spark. Quantal distribution of release flux was shown to be compatible with smaller clusters as well (Zahradnikova et al., 2010a); however, quantitative correspondence with experiments has not been tested.
In this work, the model of calcium spark (Zahradnikova et al., 2010b) has been extended and used to explore the consequences of variability of RyR clusters and their arrangement into CRUs and calcium release sites that emerged from the recent ultrastructural studies (Baddeley et al., 2009;Hayashi et al., 2009).

CONSTRUCTION OF RyR CLUSTERS, VIRTUAL CRUs, AND VIRTUAL CALCIUM RELEASE SITES
The size of the simulated RyR clusters was determined, as previously described (Zahradnikova et al., 2010a), by drawing random numbers from the bi-exponential size distribution (Baddeley et al., 2009) with the following parameters: where n RyR is the number of RyR channels in the cluster (cluster size), and p(n RyR ) is the probability that a randomly chosen cluster will have a size n RyR . A predefined number (1, 2, 3, 5, or 10) of individual RyR clusters was randomly combined into virtual calcium release units (vCRUs), which were considered to share a common dyadic space, i.e., the RyRs in a single vCRU sensed the same cytosolic and luminal environment (Figure 1Ab). One such vCRU was assumed to form a virtual calcium release site (vCRS). Alternatively, a predefined number of vCRUs (1, 2, 3, 5, or 10) consisting of a single cluster were combined into a vCRS. In this case, individual vCRUs constituting the vCRS were considered to behave independently, i.e., to have separate dyadic spaces (Figure 1Ac).

FREQUENCY OF SPONTANEOUS CALCIUM SPARKS
The frequency of spontaneous calcium sparks, F spark , was estimated according to Zahradnikova et al. (2010b):  (Zahradnikova et al., 2010b), in which the transitions between states differing in the occupation of the activation sites are described by the scheme in Figure 1B, and the channel is additionally regulated by an independent Mg 2+ inhibition site. Thus, the open probability of the RyR is: [Mg] -free Mg 2+ concentration.
where [Ca] and [Mg] are the free cytosolic concentrations of Ca 2+ and Mg 2+ , and the equilibrium constants K x and allosteric factors f x are described and the values of the parameters are given in Table 2.
Calcium spark frequency of a set of random RyR clusters was simulated using the formula derived from Eq. 2 for non-uniform cluster sizes, where N is the number of simulated vCRSs (n RyR ) i is the number of RyRs in the i th simulated vCRS and the other symbols were explained previously. www.frontiersin.org

THE NUMBER OF OPEN RyRs
The distribution of the number of open RyRs in a virtual CRU was calculated using the binomial equation (Zahradnikova et al., 2010b): where n RyR is the number of RyRs in the given virtual CRU, n O is the number of activated/open RyRs in the vCRU, and P(n O ) is the probability that n O out of n RyR RyRs will be open during the spark.
The values of p A that optimized the distribution of the number of RyRs for the whole set of simulated CRUs were determined by minimizing the sum of squares of the differences between the predicted probabilities and those observed by Wang et al. (2004): where i is the number of open RyRs, P(i) was determined by applying Eq. 4 to all vCRUs, P i is the probability of occurrence of the ith quantal level, i.e., the number of open RyRs, determined by Wang et al. (2004), and n is the number of RyRs in the largest of the simulated vCRUs. The goodness of fit for each simulation was evaluated by the χ 2 test (Press et al., 1992), using the method of Landau and Páez (1997) to determine experimental variance, as described previously (Zahradnikova et al., 1999).

PARAMETERS OF Mg 2+ UNBINDING
The parameters of Mg 2+ unbinding from the activation site of the RyR channel were determined from activation probability, p A , as described in Zahradnikova et al. (2010b). Briefly, for each simulation, the dependence of the probability of Mg 2+ dissociation (p dis ) during the spark on p A was calculated from the steady-state open probabilities of all ionic forms of RyR, defined by the number of bound Mg 2+ and Ca 2+ ions, using the parameters of the aHTG model determined for the respective simulation. The value of p dis together with an estimate of the time to the peak release flux (TTP) then yielded an estimate of Mg 2+ dissociation rate, k Mg dis , as described previously (Zahradnikova et al., 2010b): where all parameters have already been explained. Derivation of equations was performed in Mathematica (Ver. 8.0, Wolfram Research). Analyses were performed in Mathematica (Ver. 8.0, Wolfram Research) and Origin (Ver. 8.0 SR6, OriginLab).

CONSTRUCTION OF THE MODEL
A model of a segment of a cardiac myocyte, containing 100 calcium release sites (Figure 1Aa) was designed as follows: A set of 1000 RyR clusters containing altogether 14095 RyRs was generated by randomly drawing 1000 numbers from a bi-exponential distribution of cluster sizes (Eq. 1). The resulting set of random "RyR clusters" was randomly recombined into two sets of vCRSs with the same size distribution. In the first set, each vCRS contained a single supercluster CRU, into which the RyR clusters were combined, as depicted in Figure 1Ab. In the second set, the vCRS contained RyR clusters of the same size, each cluster forming a separate single-cluster CRU (Figure 1Ac).
The size distribution of the vCRSs containing 1, 2, 3, 5, or 10 clusters is shown in Figure 2 and quantitatively characterized in Table 3. For vCRSs composed of a single cluster, the size distribution was in a good agreement with the experimental data (Baddeley et al., 2009). In multi-cluster vCRSs, the minimum number of RyRs ranged from 1 to 51 and their maximum number ranged from 180 to 271, respectively, while the mean size ranged from 14 to 141 RyRs and the median size ranged from 7 to 136.5 RyRs ( Table 3). With increasing number of clusters per vCRS, the size distribution of the vCRSs gradually changed from bi-exponential to nearly Gaussian. With increasing number of clusters per vCRS, the total number of vCRSs decreased (from 1000 to 100), which was reflected in the increased variance in the size distribution histograms. When the total number of the generated vCRSs was more than 100, i.e., for n clust < 10, a set of 100 vCRSs was randomly selected for further simulation.
In all models it was assumed that all RyRs in a CRU are equivalent in the sense that they share the cytosolic as well as the luminal environment, and therefore they have the same steadystate open probability, p O , and the same probability of being activated during a spark, p A . In other words, the cytosolic and luminal environments of the RyR channels in a single CRU were considered as ideally mixed compartments. The vCRSs with multiple CRUs (Figure 1Ac) emulate the effect of insufficient spatial resolution of confocal microscopes. The probability of RyR activation in any CRU was not affected by the activation of RyRs in the remaining CRUs. Individual vCRSs were fully independent as well.

CALCIUM DEPENDENCE OF SPARK FREQUENCY
The model of calcium spark frequency is based on the idea of triggering of a spark by the first spontaneous opening of single RyR channel in the calcium release site, arising from the resting open probability of RyRs (Zahradnikova et al., 2010b). In this model, all RyR channels in the vCRS are under resting conditions before activation of a calcium spark. The activation of a calcium spark corresponds to the stochastic opening of the first RyR in the vCRS. As a result, the frequency of sparks depends only on the total number of RyRs in the vCRS (Eq. 2) and not on the number and composition of vCRUs within the vCRS.
To adjust the constructed segments to real myocytes, we tuned their calcium dependence of calcium spark frequency to that measured by Lukyanenko and Gyorke (1999). The experimental calcium dependence of spark frequency (open symbols, Figure 3) was approximated by Eq. 2, into which Eq. 3 was substituted for P O . To obtain the best approximation of the data by the models, the parameters describing interaction of Mg 2+ with the RyR activation sites (i.e., K Mg and f Mg ) and the refractory period (t rp ) were optimized. A good agreement between the theoretical calcium dependence (Eq. 2) and experimental data could be obtained

FIGURE 2 | Simulation of RyR distribution into clusters and virtual calcium release units (vCRUs).
Distribution of the size of release units (in number of RyRs per release unit, n RyR ) composed of (top to bottom) 1, 2, 3, 5, and 10 random clusters. The number of RyRs in clusters was distributed bi-exponentially according to Baddeley et al. (2009). with the mean (n RyR ) as well as with the median (ñ RyR ) number of RyRs per vCRS (albeit with slightly different optimized values of K Mg , f Mg , and t rp ). The remaining gating parameters were kept at their original values, estimated by an independent method (Zahradnik et al., 2005). The adequacy of the optimized parameters K Mg , f Mg and t rp was then tested by comparing spark frequency in the set of vCRSs, calculated using Eq. 4, with the experimental data. It turned out that to obtain a good agreement between the simulation with each set of vCRSs and the experimental data, the median (small symbols in Figure 3) but not the average of the number of RyRs per vCRS had to be used as the descriptor of vCRS size in Eq. 2. The use of the mean number of RyRs per vCRS led to disagreement between the experimentally measured and the simulated calcium spark frequencies (not shown), most probably due to the skewness of the distribution of vCRS sizes. This was especially apparent at small vCRS sizes, in which there was a large difference betweenn RyR andñ RyR ( Table 3). The median number of RyRs was an excellent predictor even for the case of the highly asymmetric bi-exponential distribution at n clust = 1. The characteristics of Mg 2+ binding to the RyR for vCRS models containing a different number of clusters (n clust ) were very sensitive to n clust (Figures 4A,B; Table 4). The value of K Mg was smaller by more than 50% when n clust was increased from 1 to 10. For n clust > 3 the optimal allosteric factor was larger than 1, as estimated previously for CRSs containing 182 RyRs (Zahradnikova et al., 2010b), implying a small negative allosteric effect of Mg 2+ on RyR open probability. For n clust ≤ 3, the value of f Mg was less than 1, implying a small positive allosteric effect of Mg 2+ on RyR open probability. The estimate of the refractory period t rp was not sensitive to changes in n clust (Figure 4C).
The Mg 2+ -binding properties of the RyR had a significant impact on RyR open probability in the whole range of cytosolic Ca 2+ , as documented in Figure 5. Here the simulated calcium dependence of open probability is shown superimposed on the published experimental data obtained in the presence of 3 mM ATP and defined concentrations of free Mg 2+ (Zahradnikova et al., 2003). The lower f Mg and the higher K Mg values at n clust of 1 and 2 led to a significantly weaker inhibition of RyR open probability than experimentally observed, especially at higher [Mg 2+ ]. The best correspondence between data and simulation occurred for n clust = 3-10. Open symbols with error bars represent data on permeabilized myocytes (Lukyanenko and Gyorke, 1999). The two indistinguishable lines represent the best fits by Eq. 2 with parameters shown in Tables 1, 2, and 4 for the optimal aHTG model for vCRSs composed of 1 and 10 random clusters with the size distribution shown in Figure 2. Small red and blue symbols represent spark frequency simulated using Eq. 4 for a set of 100 random vCRUs composed of 1 and 10 clusters.

DISTRIBUTION OF CALCIUM RELEASE FLUX
The distribution of spark calcium release flux was calculated for all simulated vCRSs using Eq. 5. When vCRSs were composed of Since under all tested conditions the resting open probability of the RyR was ∼1 × 10 −4 , the probability that a resting RyR from another cluster of the same vCRS will open during the spark was less than 0.0035, which would result in a 0.35% excess of RyR openings, less than the estimated error of approximation of the data (Wang et al., 2004) by the best model (0.95%). Therefore the contribution of resting RyR openings to the distributions of calcium release flux was not taken into account. The activation probability p A for each set of vCRSs with a given n clust was estimated by approximating the distribution of n O (Eq. 5) to the experimental data published by Wang et al. (2004). The optimal distributions are shown in Figure 6 together with the experimental distribution. The values of p A providing the optimal fit and the resulting average numbers of open RyRs per spark are given in Table 5. For vCRSs composed of a single multi-cluster vCRU (Figure 6A), it is apparent that for low n clust of 1 or 2 the theoretical distributions did not describe the data adequately (χ 2 test, p < 0.001). A satisfactory description of the data was obtained for n clust = 3-10 (χ 2 test, p > 0.1). If the vCRS was composed of a set of single-cluster, non-interacting vCRUs, the theoretical distributions of the spark calcium release flux did not depend on n clust ( Figure 6B) and were not adequate (χ 2 test, p < 0.001).

THE RATE OF Mg 2+ DISSOCIATION FROM THE ACTIVATION SITES
If the clusters in the vCRS formed a single vCRU, the probability of recruitment of additional RyRs after activation of the first RyR  (p A ) decreased with increasing n clust . We tested the extent to which these changes were due to the changes in the dissociation of Mg 2+ from the activation sites. The estimated values of p A ( Figure 7A) were converted into values of probability of Mg 2+ dissociation from the activation site during the spark, p dis (Figure 7B), using the procedure described in the Supplemental data of Zahradnikova et al. (2010b). These were further converted to Mg 2+ dissociation rates (Figure 7C), assuming that Mg 2+ dissociation has to occur not later than at the time to the peak of the spark (6.67 ms, Eq. 6).
The values of all three parameters decreased with n clust (circles); the models that did not correspond to the observed distributions of calcium release flux ( Figure 6A) provided considerably higher values (open symbols) than the accepted models (crossed symbols). If the clusters in the vCRS formed separate vCRUs, the probability of recruitment of additional RyRs into sparks after activation of the first RyR did not change considerably with n clust (squares in Figure 7). However, neither of these models corresponded to the experimental distributions of calcium release flux ( Figure 6B). Additionally, variation of the number of clusters present in one vCRS requested changes in the calcium dependence of RyR activation to follow the calcium dependence of RyR open probability at different Mg 2+ concentrations ( Figure 5). Specifically, with increasing number of clusters per vCRS, the RyRs had to undergo stronger inhibition by Mg 2+ at their activation sites (smaller K Mg and larger f Mg ). The extent, to which the changes in K Mg , f Mg , and p A are due to changes of k Mg dis , is examined in Figures 7D-F

DISCUSSION
This study reveals that not only the gating of RyR channels, but also their organization into release sites plays a role in the cardiac mechanism of excitation-contraction coupling. This would be difficult to reveal experimentally. www.frontiersin.org The relationships between the organization of RyRs into calcium release sites and the properties of calcium sparks were explored using extensions of our previously described models of calcium spark and RyR gating (see Materials and Methods). The results have shown that for a range of calcium release site sizes (n RyR = 42.3 − 141;ñ RyR = 36 − 136.5) and size distributions, a model of calcium spark formation based on independent RyRs provides a quantitative description of the calcium dependence of calcium spark frequency (Figure 3). The two alternative hypotheses, that the confocally resolved volume comprises either a single vCRU containing three or more clusters, or several small vCRUs containing one cluster each, provided identical results, since all RyRs present in the confocally resolved volume are indistinguishable under resting conditions. Moreover, the quantal distribution of spark calcium release flux in these CRUs is in a good agreement with the experimental data (Figure 6), if the RyR clusters in the calcium release sites are organized into a single CRU, that is, all RyRs of the CRS experience the same conditions modulating their steady-state open probability and probability of activation.

THE RELATION BETWEEN CRS DISTRIBUTION, CALCIUM SPARK FREQUENCY AND RyR OPEN PROBABILITY
The parameters of Mg 2+ interaction with RyR channels (K Mg and f Mg ) were optimized for different sets of vCRS sizes. A satisfactory fit of calcium spark frequency could be obtained for all vCRSs at the expense of substantial changes in the parameters K Mg and f Mg ( Table 4).
The Ca 2+ dependence of open probability of RyRs that provided the correct calcium dependence of calcium spark frequency was compared with experimental data at different Mg 2+ concentrations ( Figure 5). The simulated inhibitory effect of Mg 2+ on the calcium dependence of RyR open probability was in better accordance with the experimental observations for larger vCRS sizes (n clust × n vCRU ≥ 3). For smaller vCRSs, the correct calcium dependence of calcium spark frequency could be obtained only if RyRs inhibition by Mg 2+ at the RyR activation sites was weaker (K Mg = 150-190 µM) than experimentally estimated (61-100 µM; Gusev and Niggli, 2008;Zahradnikova et al., 2010b;Tencerova et al., 2012) except for Zahradnikova et al. (2003; 250 ± 150 µM).

QUANTAL DISTRIBUTION OF CALCIUM RELEASE FLUX
Unequivocal results regarding the composition of vCRSs were provided by analysis of the quantal distribution of calcium release flux in evoked calcium sparks. The theoretical predictions of the models corresponded with experimental data (Wang et al., 2004) only for vCRSs consisting of a single larger vCRU (n clust > 3). For vCRSs consisting of a single smaller vCRU (n clust = 1, 2), and for vCRSs composed of up to ten small vCRUs (n clust = 1), there was an excess of predicted sparks with high calcium release flux (n O > 7) and a deficit of the predicted sparks with average calcium release flux (n O = 2-3), which resulted in unacceptably high deviation of the model from the data. Since independent recruitment of RyRs after commencement of the spark, indicated by the binomial distribution of their calcium release flux, requires their equal access to the calcium trigger, all clusters in the vCRS had to communicate via the shared cytosolic/luminal environment, i.e., to consist of a single multi-cluster vCRU.

THE RATES OF Mg 2+ DISSOCIATION FROM THE RYR
The rate of Mg 2+ dissociation obtained by optimizing the quantal distribution of calcium release flux was strongly dependent on the assumed vCRS size, and ranged from 48 to 224 s −1 .A satisfactory accordance with experimental data was obtained only for vCRS size of 3-10 clusters, yielding k Mg dis = 48 − 74 s −1 . These values are comparable with or higher than the rates of dissociation of Mg 2+ from the Mg 2+ -parvalbumin complex (4-33 s -1 ; Permyakov et al., 1987) or from its complex with troponin-C or with the C-terminal fragment of troponin-C (5 and 0.7 s −1 ; Rosenfeld and Taylor, 1985). At the same time, they are somewhat lower than the rates of Mg 2+ dissociation from ATP (150 s −1 ; Baylor and Hollingworth, 1998). A relatively low value of k Mg dis may also be inferred from the strong effect of low Mg 2+ concentrations on activation of RyR channels by brief Ca 2+ stimuli (Zahradnikova et al., 2003). Therefore, the values of k Mg dis estimated for n clust of 3-10 may be considered plausible.

CONCLUSION
The modeling of the calcium release sites showed that RyRs with gating characteristics conforming to their behavior in lipid bilayers, and organized in clusters with size distribution revealed by experiments, provide a very satisfactory description of calcium spark properties in both, spontaneous and triggered sparks, collected in independent experiments. According to the best models, RyR clusters associate into CRUs composed of at least 3 clusters and containing at least 40 RyRs on average, which share the same dyadic space and environmental conditions. Moreover, these results reconcile the observed small cluster size and somewhat irregular placement of RyRs in the release sites with the previous estimates of the amount of RyR protein (Bers and Stiffel, 1993), volume density of calcium release sites (Chen-Izu et al., 2006), and the extent of calcium release in rat cardiac myocytes (Zahradnikova et al., 2010b).