On the Adjacency Matrix of RyR2 Cluster Structures

In the heart, electrical stimulation of cardiac myocytes increases the open probability of sarcolemmal voltage-sensitive Ca2+ channels and flux of Ca2+ into the cells. This increases Ca2+ binding to ligand-gated channels known as ryanodine receptors (RyR2). Their openings cause cell-wide release of Ca2+, which in turn causes muscle contraction and the generation of the mechanical force required to pump blood. In resting myocytes, RyR2s can also open spontaneously giving rise to spatially-confined Ca2+ release events known as “sparks.” RyR2s are organized in a lattice to form clusters in the junctional sarcoplasmic reticulum membrane. Our recent work has shown that the spatial arrangement of RyR2s within clusters strongly influences the frequency of Ca2+ sparks. We showed that the probability of a Ca2+ spark occurring when a single RyR2 in the cluster opens spontaneously can be predicted from the precise spatial arrangements of the RyR2s. Thus, “function” follows from “structure.” This probability is related to the maximum eigenvalue (λ 1) of the adjacency matrix of the RyR2 cluster lattice. In this work, we develop a theoretical framework for understanding this relationship. We present a stochastic contact network model of the Ca2+ spark initiation process. We show that λ 1 determines a stability threshold for the formation of Ca2+ sparks in terms of the RyR2 gating transition rates. We recapitulate these results by applying the model to realistic RyR2 cluster structures informed by super-resolution stimulated emission depletion (STED) microscopy. Eigendecomposition of the linearized mean-field contact network model reveals functional subdomains within RyR2 clusters with distinct sensitivities to Ca2+. This work provides novel perspectives on the cardiac Ca2+ release process and a general method for inferring the functional properties of transmembrane receptor clusters from their structure.


Introduction
Mechanical contraction of the heart occurs as a result of intracellular Ca 2+ release in cardiac myocytes. L-type Ca 2+ channels (LCCs) and a packed cluster of up to 100 Ca 2+ -sensitive Ca 2+release channels [1,2], known as ryanodine receptors (RyR2s), are co-located at discrete subcellular junctions within the cell (Fig 1A). These Ca 2+ release units (CRUs) are formed by deep invaginations of the cell membrane containing LCCs, known as transverse-tubules (TTs), and the junctional sarcoplasmic reticulum (JSR) membrane, a cisternal sheet containing the RyR2s that wraps around the TT to form a narrow subspace * 15 nm in width. During excitationcontraction coupling (ECC), electrical stimulation increases the probability of LCC openings and influx of Ca 2+ into these subspaces. Binding of Ca 2+ to the closely-apposed RyR2s [3] increases their open probability and release of Ca 2+ from JSR stores in a process known as Ca 2+ -induced Ca 2+ release (CICR). Further Ca 2+ release from RyR2s activates surrounding RyR2s via a local rise in subspace Ca 2+ concentration. Understanding this process is critical to our understanding of cardiac physiology in health and disease.
Isolated release events known as Ca 2+ "sparks" underlie the cell-wide release of Ca 2+ that occurs on every heartbeat when the RyR2s are activated by the opening of voltage-sensitive Ca 2+ channels. Ca 2+ sparks are also observed in resting myocytes when initiated by the spontaneous opening of a single RyR2 that then probabilistically triggers Ca 2+ release from the rest of the cluster. Here, the probability that an RyR2 is in an open state at any point in time will be referred to as its open probability, and the probability that a sufficiently large percentage of the RyR2s open for a spark to occur will be referred to as spark probability (p S ).
Spark probability is an important physiological parameter that in part controls the frequency of sparks [4,5]. There is significant experimental evidence that not all junctional RyR2 openings result in Ca 2+ sparks [6][7][8][9]. These non-spark openings may in part be attributed to non-junctional RyR2s located outside the release site [10,11]. However, mathematical modeling suggests that junctional RyR2 openings can fail to trigger Ca 2+ sparks and are sufficient to account for the non-spark openings [4,5,12,13]. In support of this, recent experiments using Ca 2+ nanosensors targeted to the release site have demonstrated locally elevated Ca 2+ concentration in the subspace due to spontaneous junctional RyR2 openings [14,15]. This is also consistent with the observation that the majority of Ca 2+ released via non-spark openings is extruded through the Na + /Ca 2+ exchanger, which is localized near the junctions [15]. Therefore there is compelling evidence suggesting that Ca 2+ spark initiation is likely a probabilistic process.
Many studies have implicated Ca 2+ sparks in heart disease. For example, spark frequency is increased in heart failure [7,16], which is associated with decreased JSR Ca 2+ content and thus impaired contractile function [17]. Ca 2+ sparks may also cause spontaneous Ca 2+ waves [18,19] that promote cellular arrhythmias [20]. Heterogeneity in the Ca 2+ sensitivity of RyR2 clusters has also been implicated in the occurrence of arrhythmic Ca 2+ alternans [21]. Therefore factors that influence Ca 2+ spark probability are likely to be involved in mechanisms driving pathological cardiac dysfunction.  49. After a single channel is opened at t = 0, it probabilistically succeeds (blue) or fails (red) to trigger a spark. Inset shows the spark initiation phase. (D) Relationship of cluster size (n, triangles) and maximum eigenvalue (λ 1 , circles) with Ca 2+ spark probability in the 3D spark model [4]. Two embedded cluster lattices are shown to illustrate that RyR2 clusters can exhibit equal Ca 2+ spark probability despite differences in the number of channels due to differences in shape.
Advancements in super-resolution imaging techniques have enabled the study of nanoscale receptor organization in a variety of cell types [22,23]. In cardiac myocytes, stimulated emission depletion (STED) microscopy has been applied to study TT remodeling in heart failure at nanometer resolution [24]. Within the cardiac Ca 2+ release sites, RyR2s are known to form tightlypacked clusters from electron microscopy studies in vivo [1] (Fig 1B), which is supported by the observation that the channels organize into packed lattices with * 31 nm spacing in vitro [25]. Super-resolution imaging studies of RyR2 clusters using super-resolution microscopy in cardiac myocytes have revealed that they are heterogeneous in size and shape [2].
We recently developed a three-dimensional, biophysically-detailed model of cardiac Ca 2+ release events, which we will refer to as the 3D spark model. Fig 1C shows example simulations, in which a single open channel either fails or succeeds to activate the remaining RyR2s. We obtained realistic RyR2 clusters obtained using STED microscopy and showed that the precise spatial arrangement of RyR2s critically influences spark probability [4]. Larger, more compact clusters exhibited higher spark probability than smaller, fragmented ones. Representing RyR2 clusters as a two-dimensional lattice, we found that the maximum eigenvalue (λ 1 ) of the lattice's adjacency matrix predicts the Ca 2+ spark probability of the cluster. Properties of the eigenvalues and eigenvectors of a graph's adjacency matrix have been widely studied, and λ 1 in particular is known to be a measure of the interconnectedness of the graph [26,27]. λ 1 was found to be a more accurate predictor of spark probability than is the total number of channels (n), which does not consider structural aspects of the cluster. Fig 1D shows the relationship between spark probability, n, and λ 1 for two RyR2 clusters analyzed previously [4]. The two embedded cluster lattices have nearly equal spark probability, despite one being much larger than the other. This is because the larger lattice contained four empty spaces in the interior. However, the λ 1 values were similar for these two clusters and therefore consistently correlated with spark probability.
Here we present an analytical model of the Ca 2+ release process and derive the relationship between λ 1 and spark probability. The model is applied to realistic RyR2 clusters obtained using stimulated emission depletion (STED) microscopy. We found through an eigendecomposition that some RyR2 clusters possess functional subdomains with distinct sensitivity to Ca 2+ . This work outlines a unique approach to understanding CICR and provides a theoretical framework for comparing the physiological function of protein clusters based solely on structural information.

Contact network model
In this section, we present results using a contact network (CN) model of RyR2 cluster activation where channels are coupled through local interactions with their neighbors. Contact network models are widely used to study the spread of disease due to contact between infected and susceptible individuals [28]. In our model, interactions instead arise from Ca 2+ -dependent activation due to local influx and diffusion of Ca 2+ , which causes neighboring channels to open. For simplicity, we assume that the local Ca 2+ concentration gradient near an open RyR2 declines rapidly enough in space such that only adjacent RyR2s interact [4,29,30]. Each channel transitions stochastically between open and closed states (Fig 2A). If an RyR2 channel i has Y i (t) neighboring RyR2 channels that are open, its opening rate is βY i (t), where β is a constant parameter. Therefore β is the RyR2 opening rate when one nearest neighbor RyR2 is open. Note that in the full biophysical 3D spark model, the RyR2 opening rate when all neighbors are closed is very small (* 9 × 10 −7 ms −1 ). Therefore we have taken this rate to be zero in this formulation. The value of β is varied in our analyses. The RyR2 closing rate, δ, is assumed to be a constant 0.5ms −1 . Derivation of the model and parameters are given in Methods and Models. The CN model is able to capture RyR2 gating dynamics during the initiation phase of Ca 2+ sparks. We used the Stochastic Simulation Algorithm of Gillespie [31] to simulate the stochastic CN model. Fig 2B shows traces of the number of open channels (N O ) during representative simulations of spark initiation in the 3D spark and CN models for a 7 × 7 lattice cluster. A single RyR2 is opened at t = 0, which then triggers openings of other channels. The CN model qualitatively reproduces channel gating behavior during the initiation of the spark. In the 3D spark model, Ca 2+ sparks occur with greater than 95% probability if a minimum of four channels open. Therefore, we define this as the minimum number for successful spark initiation in both models. We also assume that each RyR2 in the cluster is equally likely to open spontaneously, and so the first open channel is chosen at random.
The advantage of developing the CN model is that we can derive analytical relationships between the dominant eigenvalue of the RyR2 lattice's adjacency matrix, λ 1 , and spark probability. We show (see Eq (9) in Methods and Models) that for a deterministic mean-field approximation of the model, RyR2 open probability decays to zero when This implies that, in the mean-field approximation, δ/β is a stability threshold for λ 1 at which RyR2 activity switches from decay to growth. While it was not immediately clear how this threshold related to the behavior of the full stochastic CN model, we expected that the model would exhibit constant spark probability when λ 1 = δ/β. That is, for a set of cluster structures each with a different value of λ 1 , the spark probability would be consistent across clusters when each cluster's opening rate was set to β = δ/λ 1 . Fig 2C shows the spark probability for a collection of 107 RyR2 clusters obtained using STED microscopy (see [4] for imaging methods). For each simulation, λ 1 was computed for the cluster and β was set to the threshold value δ/λ 1 . The range of spark probabilities across all clusters was narrow (0.14±0.0078). This was also observed when using sub-threshold values β = 0.5δ/λ 1 (0.029±0.0024) and a supra-threshold values β = 2δ/λ 1 (0.28±0.012). Therefore spark probability is constant when β is scaled inversely with λ 1 . For comparison, we also plotted spark probability when β is set to a single value across all clusters ( Fig 2D). In this case, spark probability increased with λ 1 in agreement with the 3D spark model (see Fig 1D). The CN model was able to accurately predict Ca 2+ spark probability for a range of cluster geometries. We estimated the spark probabilities for a collection of 15 RyR2 clusters obtained using STED microscopy. The value of β was adjusted until the spark probabilities in the CN model correlated with those of the 3D spark model (Fig 2E). Maximal correlation was achieved for β = 0.115 (R 2 = 0.939), which gives the value of δ/β = 4.35. Note that the theoretical value of λ 1 for any cluster is bounded above by the maximum number of channel neighbors (4) [26]. Consequently, δ/β = 4.35 > λ 1 implies that the system is always sub-threshold for any cluster structure under normal physiological conditions. The CN model also predicts spark probability for different opening rates. To show this, we first estimated p S in the 3D spark model for a 7 × 7 cluster with the opening rate scaled by a constant factor. We then scaled β = 0.115 by the same factor and determined p S in the CN model. This was repeated for a range of scaling factors. Noting that the closing rates δ are the same in both models, we could directly compare p S in the two models by plotting it as a function of δ/β, where β is the scaled value. For the 3D spark model, β is the value used in the corresponding CN simulation. In both models, spark probability fell rapidly as δ/β approached λ 1 from the left before decreasing gradually to the right of λ 1 . This suggests that spark probability is more sensitive to RyR2 gating kinetics when the opening rate is elevated. From the data in this section, we concluded that the CN model is able to accurately predict p S over a range of opening rates and cluster geometries.

Ca 2+ diffusion in the CN model
Cardiac Ca 2+ release is actively regulated under normal conditions and modulated in various diseases. To study this regulation, we expanded the CN model by deriving a simple model of Ca 2+ diffusion between RyR2 Ca 2+ sources. The parameter β was estimated using this diffusion model and a model of RyR2 gating. All parameters were taken from Walker et al. [4], except for the effective Ca 2+ diffusion coefficient (d C ), which was adjusted to give β = 0.115 as determined in the previous section.
A number of signaling molecules regulate RyR2 channels, affecting their opening rate. This includes RyR2 phosphorylation by Ca 2+ /calmodulin-dependent protein kinase II (CaMKII) and protein kinase A (PKA) [32,33] and JSR Ca 2+ concentration [34]. Channel gating can also be altered under oxidative stress [35] and by genetic mutations [36,37]. As shown in Fig 3A, δ/ β is inversely proportional to the channel opening rate constant (k + ), reflecting the increased Ca 2+ spark frequency observed under such conditions [7,38,39]. Note that the closing rate is δ and therefore scales δ/β linearly. Increasing the unitary channel current (i RyR ) resulted in a decrease in δ/β ( Fig 3B). This behavior is consistent with experimental evidence [40], in which decreased i RyR resulted in lowered spark frequency.  The CN model was also sensitive to parameters affecting the diffusion of Ca 2+ ions in the release site subspace. Fig 3C shows the dependence of δ/β on d C . As d C increases, Ca 2+ ions are more likely to escape the nanodomain around the open channel, thus decreasing spark probability. Uniformly increasing the distance between the open channel pore and neighboring Ca 2+ binding site increased δ/β so as to decrease spark probability (Fig 3D).
In this section, we have used a simple diffusion model to probe the effects of perturbations to biophysical properties of the release site including the opening rate, unitary channel current, Ca 2+ diffusion coefficient, and inter-channel spacing. The CN model suggests that minor modifications to these parameters can alter the stability of the system, thus leading to significant changes in spark probability.

Linear mean-field CN model
Up to this point, we have considered spark probability when each RyR2 is equally likely to open first. An emergent property of the 3D spark model was that the probability of a spark occurring varied with the choice of initiating RyR2 [4]. Channels closer to the epicenter of the cluster were more likely to trigger sparks because they have more possible combinations of first, second, third, etc. neighbors along which channel openings could propagate. Likewise, channels on the periphery of the cluster were less likely to trigger sparks.
We To further establish the relationship between E[N O ] and p S , we compared E[N O ] to λ 1 for a broader collection of 107 clusters obtained from STED microscopy ( Fig 4B). A strong correlation between these variables was present, in contrast to the number of channels in the cluster, which was not consistently correlated with E[N O ]. Recall that these conclusions were also drawn from the data of Fig 1D for a smaller collection of clusters. Taken together, these data suggest that λ 1 and E[N O ] are both accurate predictors of p S , while by itself the number of channels without regard to relative channel locations is not.
The LCN model was used to compute the vector whose elements are the expected number of open channels given each possible initiating RyR2. We denote this vector E n O ½ (see Eq (14) (15) and (14)). RyR2. ðE n O ½ Þ j decayed more rapidly for channels j near the edge compared to those near the center, consistent with the lower peripheral spark probabilities estimated using the 3D spark model (Fig 4D). This is because channels near the edge have fewer adjacent channels to trigger, and therefore it is less likely that a second channel will open before the first closes.
Furthermore, the peripheral channels tend to have fewer second, third, etc. neighbors that can potentially be activated compared to central channels. We conclude that both the expected number of open channels in the LCN model is strongly correlated with spark probability. This fact will be used to further analyze spatial gradients in spark probability that depend on which RyR2 is opened initially.

Characterization of functional subdomains by the eigenmodes
The LCN model can be decomposed into a set of independent eigenmodes by taking the similarity transform of the adjacency matrix: A = V D V T , where V is the modal matrix whose columns are formed by a set of orthonormal eigenvectors fv 1 ; v 2 ; :::; v n g and D is a diagonal matrix of eigenvalues {λ 1 , λ 2 , . . ., λ n } in descending order. Note that A is symmetric such that V −1 = V T . The i th eigenmode is defined by the pair l i À v i , in which λ i determines the rate of decay of the eigenmode in time and the values ðv i Þ j determine the membership of channel j in the eigenmode. We derived an expression for E n O ½ as a weighted sum of the eigenvectors where are essentially equal to weighted sums of the eigenmodes. The derivation of these equations can be found in Methods and Models.
In the previous section, we presented further evidence of the relationship between λ 1 and spark probability as well as intra-cluster spatial gradients in spark probability. A natural question to then ask is: does the dominant eigenvector (v 1 ) corresponding to λ 1 give us information about these gradients? Furthermore, how significant are other eigenmodes?
The spatial distribution of E n O ½ is shown for collection of 10 RyR2 clusters in Fig 5A. We further defined c i ¼ c i ðtÞ= P j c j ðtÞ att ¼ 8 ms, which gives the fractional contribution of the i th eigenmode to N O ðtÞ. We decomposed these clusters into their eigenmodes and plotted the values of c i corresponding to the 8 greatest eigenvalues in Fig 5A. In most cases (clusters (1)-(4), (6), (7), (9)), c 1 was the only large value, implying that the dominant eigenmode characterized the behavior of the LCN model. Clusters (5), (8), and (10), however, exhibited another significant c i corresponding to a subdominant eigenmode.  Examining the dominant eigenmode of each cluster, we found that for clusters characterized by only the dominant eigenmode, there was a single locus of elevated membership in the dominant eigenvector corresponding to the channels j with greatest values E n O ½ ð Þ j (Fig 5C). Furthermore, the eigenmode's spatial gradients resembled the full solution in Fig 5A. For clusters with a subdominant eigenmode ((5), (8), (10)), however, the dominant eigenvector did not fully characterize the spatial gradients in E n O ½ . For these clusters, the subdominant eigenmode accounted for areas of high E n O ½ that were not included in the dominant eigenmode (Fig 5D). In addition, the subdominant eigenmodes were insignificant in the other clusters.
To quantitatively assess how well the dominant and subdominant eigenmodes characterize spark probability, we computed the correlation coefficients between the E n O ½ and the dominant (k = 1), the sum of dominant and subdominant (k = 2), and sum of the dominant, subdominant, and a tertiary eigenmode corresponding to the third largest c i (k = 3) (Fig 5E). Clusters well-described by the dominant eigenmode generally yielded high ρ 1 > 0.84, indicating that v 1 was sufficient to characterize the spatial gradients in spark probability. For the clusters with significant subdominant eigenmodes, ρ 1 was lower (< 0.68), and the second eigenmode was required to establish a correlation. Note that inclusion of the tertiary eigenmode did not greatly improve the correlation, suggesting that the first two eigenmodes were most significant.
In this section, we characterized the intra-cluster spatial gradients in spark probability in terms of the eigenvalues and eigenvectors of the adjacency matrix. In the majority of cases, the dominant eigenmode l 1 À v 1 was sufficient to approximate the gradients. Clusters (5), (8), and (10) of Fig 5, however, possessed secondary subdomains of channels separated from the dominant subdomain by a bottleneck (i.e. dumbbell-like morphology). These functional subdomains generally contained channels with lower spark probability than the dominant subdomain. This is consistent with Eq (2), which indicates that these secondary subdomains are also characterized by a decay rate 1/λ s > 1/λ 1 and therefore would be expected to have lower spark probability.

Perturbation analysis
It is not clear how one can determine whether a cluster is characterized by a single eigenmode or dominant-subdominant pair of eigenmodes without performing the eigendecomposition computations. For example, comparing clusters (6) and (8) in Fig 5, it is not immediately obvious why (8) requires both modes and (6) does not. To better understand this relationship, we progressively severed the connection between two functional subdomains at the top and bottom of cluster (6). In Fig 6A, we removed channels from this cluster proceeding left to right along the row of channels indicated by the dashed line in the baseline cluster. A subdominant eigenmode emerged as the channels were removed. The dominant eigenmode remained in the lower subdomain, while the subdominant eigenmode formed in the upper region. Note the formation of two disjoint subclusters in cluster (A4), which have eigenmodes similar to when connected by a single channel in (A3). The formation of a secondary subdomain is further demonstrated by an increase in the value of c i for the subdominant eigenmode (Fig 6B). In this example, the subdominant eigenmode appeared after removing only one channel and gradually became more prominent with the removal of additional channels. Therefore, the formation of a subdominant eigenmode can be quite responsive to reductions in the region dividing two possible subdomains, each distinguished by different propensities for sparks.
We next investigated how sensitive spark probability is to small changes in lattice shape. Fig  6C shows a series of clusters in which only a single channel was removed from the original cluster. As expected, removing a channel along the upper edge as in cluster (C1) where spark probability is low resulted in a small change in λ 1 (Δλ 1 ). Discarding a central channel as in cluster (C4) resulted in the greatest change. One may expect that removing channels with the higher spark probability would cause a greater decrease in λ 1 . However, the channel removed in cluster (C2) corresponded to a greater element in E n O ½ than the channel in (C3), yet the change in λ 1 was less. This is because the channel in cluster (C3) had greater membership in the dominant eigenmode. To illustrate this, we systematically removed each channel one at a time from the baseline cluster and calculated Δλ 1 . Fig 6D shows that there was a consistent relationship between Δλ 1 and the discarded channel j's corresponding element of the dominant eigenvector ðv 1 Þ j in the original cluster but not ðE n O ½ Þ j . Therefore element j of the dominant eigenvector determines the extent by which spark probability decreases when a single channel j is removed.

Discussion
We have shown in previous work that the precise structure of RyR2 channel clusters influences properties of Ca 2+ release [4]. In particular, the probability of a Ca 2+ spark occurring when an RyR2 opens spontaneously depends strongly on the arrangement of the RyR2s in the subspace. This has implications for Ca 2+ cycling in the cell, as Ca 2+ spark probability controls the  On the Adjacency Matrix of RyR2 Cluster Structures frequency of Ca 2+ sparks and the excitability of the cluster [5]. An emergent property of this biophysically-detailed model was that λ 1 is a strong predictor of Ca 2+ spark probability.
Here we presented a model similar to those used to study the spread of contagion, such as in disease epidemics [41]. In this model, a single RyR2 is opened initially, which increases the open probability of its neighbors via a local rise in Ca 2+ concentration. After deriving a linearized mean-field formulation of the system, we showed that the open probability of the RyR2s is extinguished when λ 1 < δ/β. Therefore λ 1 governs a stability threshold for spark generation. In the stochastic model, spark probability was constant across all clusters when λ 1 = δ/β. Therefore, if one compares any two different RyR2 clusters, the cluster with lower λ 1 value would need a lower δ/β ratio (i.e. greater RyR2 mean open time or opening rate) to achieve the same spark probability as the other cluster. This explains why λ 1 is correlated with Ca 2+ spark probability. Cator and Van Mieghem derived a second-order CN model, with which they showed that the true threshold for the system to exhibit exponentially long transients is in fact bounded from above by λ 1 [42]. Nevertheless, the first-order model presented here was sufficient to account for the relationship between λ 1 and spark probability.
It is known that the maximum eigenvalue of a graph's adjacency matrix is related to the number of walks on the graph [26]. Specifically, if W k is the number of possible walks of length k on a graph with n vertices, then W k % nl k 1 when k is large. Furthermore, W k is proportional to ðu T v 1 Þ 2 when k is large. Recall that this term also appears in the expression for c 1 . It is no coincidence that W k and E[N O ] are related. Intuitively, a greater number of walks implies that there are more possible contiguous sets of RyR2s along which channel openings can propagate. This is essentially the underlying relationship between λ 1 and Ca 2+ spark probability. An eigendecomposition of the CN model further identified RyR2 subdomains characterized by different spark probabilities, as observed in the 3D spark model. Secondary subdomains with lower spark probability were found in clusters containing two groups of channels separated by central narrow regions * 2 -3 channels in width. This lends meaning to the eigenvectors of the model, which define the membership of the RyR2s to each functional subdomain. Interestingly, v 1 is a known measure of vertex centrality [43], which means that the proportion of all possible walks of length k beginning at vertex j is ðv 1 Þ j =ðu T v 1 Þ when k is sufficiently large. This implies that the elements of v 1 indicate the relative number of lattice walks beginning at each channel. Our results suggest that this is approximately true for clusters characterized by the dominant eigenmode, as channels with greater values of ðv 1 Þ j had higher spark probability. Because we consider the transient behavior of channel gating during a fixed time window, the assumption that k is large (i.e. t is large) may not hold, thus explaining why a subdominant eigenmode was observed for some clusters.
Our results indicate that the system is near the threshold under normal conditions, as the ratio δ/β is close to the threshold λ 1 . Therefore, small changes to β can greatly change the qualitative behavior of the system. Using a simple Ca 2+ diffusion model, we determined that spark probability is sensitive to changes in biophysical parameters.
RyR2 open probability is modulated by a variety of factors including phosphorylation [32,33], JSR Ca 2+ concentration [34], oxidative stress [35], and genetic mutations [36,37]. Most of these increase the opening rate of the channel and cause elevated Ca 2+ spark frequency. Our recent work [4] and others [30,44] have shown that RyR2 regulation by JSR Ca 2+ concentration is not necessary for spark termination. Rather, depletion of the JSR Ca 2+ stores causes a sufficient decrease in unitary RyR2 current such that the channel openings are not sustained. This mechanism is supported in the present model as well, as shown by the sharp increase in δ/ β as i RyR is decreased (see Fig 3B), as it would be due to reduction of the Ca 2+ concentration gradient from inside the JSR to the subspace when RyR2s open.
A recent imaging study by Asghari et al. [45] observed regulation of RyR2 cluster structure. The authors reported RyR2s clusters in dense side-by-side lattices, as assumed in the present study, as well as checkerboard-like arrangements with greater spacing of * 37 nm compared to the baseline of 31 nm. Increasing channel spacing uniformly caused an increase in δ/β to 6.3 at 37 nm (see Fig 3D). Note that for any graph whose vertices have a maximum of m neighbors, λ 1 < m [27]. Therefore λ 1 < 4 for cluster lattices. This suggests that any cluster in the checkerboard arrangement would be unlikely to exhibit Ca 2+ sparks in the absence of other changes. Interestingly, checkerboard spacing was observed upon channel phosphorylation or after decreasing the cytosolic Mg 2+ concentration, both of which increase RyR2 open probability [33,34]. Therefore the increase in inter-channel spacing may counteract the effects of these conditions.
We maintained focus on the relevance of cluster morphology to Ca 2+ spark probability when a single RyR2 opens spontaneously. Ca 2+ release can also be triggered following electrical excitation of the cell due to Ca 2+ influx through apposing LCCs located on the transverse tubule. Note modeling studies suggest that coupling fidelity between LCCs and RyR2s is still strong despite low spark probability [4,5,13]. This is because although LCC mean open time is shorter (* 0.5 ms), unitary LCC current is approximately the same as the RyR2, there are usually several LCCs per RyR2 cluster (the ratio of RyR2s to LCCs is 4-10 [46]), and LCC openings are synchronized upon membrane depolarization to drive local buildup of Ca 2+ .
The study of sub-cellular structure using super-resolution techniques requires careful interpretation of the raw image data. In this study, we generated RyR2 cluster lattices based on fluorescence intensity using a uniform thresholding algorithm. Intensities at or above the 95 th percentile were interpreted to represent the RyR2 positions over the entire image. To assess uncertainty in the results with respect to our choice of threshold, we analyzed a single set of STED images using both the 95 th and 98 th percentile thresholds. At the higher threshold, more of the fluorescence signal is filtered out and thus the clusters contain fewer RyR2s. This resulted in lower values of λ 1 and decreased p S (Fig 7A). The large differences in spark probability after using the higher threshold highlight the sensitivity of the model to the image processing methods. Nevertheless, there was still a strong correlation between p S and λ 1 when using the higher threshold ( Fig 7B). Consequently quantitative prediction of spark probability applying λ 1 requires consistent interpretation of super-resolution imaging data and in addition benefits from an incremental alteration of image analysis parameters if possible.
We did not consider weaker interactions between RyR2s such as those between diagonallyadjacent neighbors. This results in an underestimation of the open probabilities. We also did not consider clusters with heterogeneous inter-channel spacing as observed in Asghari et al. [45]. We also only considered single connected clusters containing no gaps that divide the cluster into separate subclusters. We assumed that the Ca 2+ concentration gradient surrounding an open RyR2 declines sufficiently rapidly such that a negligible Ca 2+ concentration is sensed in nearby subclusters, and therefore spark initiation occurs independently. These limitations could be overcome by using a distance matrix or diffusion model as in [47] to compute inter-RyR2 Ca 2+ coupling. In addition, the LCN model is known to deviate most from the exact model near the stability threshold (δ/β % λ 1 ) [48]. Note it has been shown that the solution of the mean-field CN model is an upper-bound on the true probabilities [49], and although higher-order approximations do exist [42] we chose the first-order mean-field approximation for its simplicity and analytical tractability.
The CN model may also be applied to similar biological systems. It may be adapted to study Ca 2+ release triggered by an LCC. The spark probability would be related to the coupling fidelity between LCCs and RyR2s. This model could be used to analyze the arrangement of LCCs as such experimental data become available. It may also be applied to future imaging studies to compare RyR2 cluster morphology to, for example, identify interspecies variability or remodeling in heart disease. For example, reduced RyR2 cluster sizes and fragmented JSR morphology have been observed in mouse models of catecholaminergic polymorphic ventricular tachycardia [50]. Inositol trisphosphate receptors (IP3Rs) located in the sarcoplasmic reticulum are known to aggregate into small clusters that exhibit similar release events known as Ca 2+ "puffs," and recent work has implicated cluster size in release extent [51] and trigger probability [52]. The present models could be used to compare IP3R cluster geometries like those reported in a recent study [53]. In skeletal muscle, Ca 2+ release is coordinated mainly by physical LCC-RyR1 and RyR1-RyR1 interactions [54]. Imaging studies have observed that RyR1 clusters in slow-twitch muscle fibres were typically smaller and more fragmented than in fasttwitch muscle [1,55]. The model presented here could be used to relate these observations to known differences in the Ca 2+ release properties of these cell types.
More generally, the model could be adapted to complement super-resolution imaging studies of a wide range of receptors that form similar supramolecular clusters in other cell types [22,23]. A general theoretical model has suggested that clusters of ligand-activated receptors behave cooperatively [56]. Examples from neurons include include synaptic microclusters of syntaxin 1 [57], acetylcholine receptor complexes at the neuromuscular postynapse [58], and On the Adjacency Matrix of RyR2 Cluster Structures rings of AMPA receptors found in spiral ganglion neurons [59]. Another example are immunoreceptors [60], which form clusters to amplify signal initiation and transduction, perhaps by decreasing the effective dissociation constants of ligands and downstream effectors [61]. Furthermore, Greenfield et al. employed super-resolution techniques to investigate the spatial organization of receptors involved in bacterial chemotaxis [62]. These receptors form clusters in the cell membrane and, similar to RyR2s, exhibit cooperative interactions that enhance sensitivity to low chemical signals.
This work presents a new perspective on cardiac calcium release and, more generally, highlights the relevance of subcellular variability in microdomains for the study of multi-scale biological systems.

Ethics statement
All animal procedures were reviewed and approved by the Institutional Animal Care and Use Committee at University Medicine Göttingen Zentrale Tierexperimentelle Einrichtung (ZTE). Animal sacrifice was applied as described in Wagner et al. [24].

Contact network model
Contact process models have been widely studied for their use in modeling disease and computer virus spread (see Keeling and Eams [41] for a review). In the present work, the CN model represents the RyR2 channel gating of a cluster of n channels. We will restrict ourselves to clusters that are connected, i.e. there are no separate islands of channels. The model is composed of a set of n random variables X i (t) = 1 if channel i is open at time t and 0 otherwise. If the channel is open, the probability that it closes within an infinitesimal time step dt is given by δdt, where δ = 0.5ms −1 is constant. If channel i is closed, it transitions into the open state in time dt with probability βY i (t)dt, where Y i (t) is the number of open adjacent channels. β is a constant given by β = k + C η , where k + = 1.107 × 10 −4 μM −η s −1 is the opening rate constant, C is the local elevation of Ca 2+ concentration caused by an open neighbor, and η = 2.1 is the Hill coefficient for Ca 2+ binding [4].
The adjacency matrix A is defined as an n × n matrix, where element (A) ij = 1 if channels i and j are adjacent, and 0 otherwise. The number of open adjacent channels is then given by Y i (t) = ∑ j (A) ij X j (t). Let p i (t) = P(X i (t) = 1), the probability that channel i is open at time t, which obeys the equation [48] ðAÞ ij X j ðtÞ À dX i ðtÞ: ð4Þ The entire system can be more compactly represented as the matrix equation where pðtÞ ¼ ½p 1 ðtÞ; :::; p n ðtÞ, u is the all-one-vector, XðtÞ ¼ ½X 1 ðtÞ; :::; X n ðtÞ, and I is the identity matrix. The system is therefore described by a set of n coupled stochastic differential equations, whose solution is analytically intractable. We simulated the CN model using the Gillespie algorithm [31]. Spark probability in the CN model was estimated by running an ensemble of 10,000 simulations per data point.
where i RyR = 0.15 pA is the unitary current of a single channel, z = 2 is the valency of Ca 2+ , F is Faraday's constant, d C is the effective diffusion coefficient of Ca 2+ in the release site subspace, and r = 31 nm is the distance between the open channel pore and neighboring Ca 2+ binding site. The diffusion coefficient for Ca 2+ in the subspace is unknown, though estimates for d C in the cytosol range from 100 to 600 μm 2 s -1 [64]. Ca 2+ buffering molecules, electrostatic interactions with the membrane, and tortuosity imposed by the large RyR2 channels can affect the motion Ca 2+ ions [65]. In light of these factors, the value of d C was adjusted from 250 to 146 μm 2 s -1 to obtain the nominal value of β = 0.115 that yields accurate spark probabilities (see Results).

Linear mean-field contact network model
A common approach is to derive a mean-field approximation of the first moment of X i (t) by assuming that the higher moments are equal to 0 [48]. This yields a set of non-linear ordinary differential equations dpðtÞ dt ¼ bdiagfu À pðtÞgA À dI pðtÞ; ð7Þ where pðtÞ is now the vector of mean-field open probabilities. This non-linear system is difficult to analyze analytically [48]. We further simplify the model by linearizing the equations about p ¼ 0 [28] dpðtÞ dt ¼ bA À dI ð Þ pðtÞ: We refer to this as the linearized mean-field CN (LCN) model, which is amenable to the tools of linear systems theory. Note that the system is stable if and only if the maximum (dominant) eigenvalue of β A − δ I, given by βλ 1 −δ, is less than 0, or where λ 1 is the maximum (dominant) eigenvalue of A. Therefore, if λ 1 < δ/β, the open probabilities in the LCN decay to 0. Otherwise, pðtÞ is unbounded as t ! 1. While physically meaningless, this result implies that the open probabilities increase when most channels are closed, or pðtÞ % 0 (near the origin of linearization).
The eigendecomposition of A is given by where V is the modal matrix with columns formed by the orthonormal eigenvectors fv 1 ; :::; v n g of A, and D is a diagonal matrix of the eigenvalues {λ 1 , . . ., λ n } in descending order.
Note that A is symmetric and therefore V −1 = V T . Combining Eqs (8) and (10) gives which can be rewritten as the summation The solution of this system is given by We refer to the eigenmodes as the eigenvalue-eigenvector pairs l i À v i . Note that pðtÞ is essentially a sum of the eigenmodes. If the initial probability distribution pð0Þ ¼ av i for some constant α, then pðtÞ / v i for all t. In other words, the trajectory of the system will be entirely characterized by the i th eigenmode. In general, the contribution of the i th eigenmode is determined by the weight v T i pð0Þ and a time-dependent exponential factor with time constant 1/(βλ i − δ).
We define E n O ðtÞ ½ as the vector whose elements E n O ðtÞ ½ ð Þ i give the expected number of open channels at time t given that channel i is open initially. This is computed by taking the sum of the elements of pðtÞ in Eq (13) We assume that in a resting RyR2 cluster, every channel experiences the same Ca 2+ concentration and therefore is equally likely to initiate a spark.