Nanoscale organization of ryanodine receptor distribution and phosphorylation pattern determines the dynamics of calcium sparks

Super-resolution imaging techniques have provided a better understanding of the relationship between the nanoscale organization and function of ryanodine receptors (RyRs) in cardiomyocytes. Recent data have indicated that this relationship is disrupted in heart failure (HF), as RyRs are dispersed into smaller and more numerous clusters. However, RyRs are also hyperphosphorylated in this condition, and this is reported to occur preferentially within the cluster centre. Thus, the combined impact of RyR relocalization and sensitization on Ca2+ spark generation in failing cardiomyocytes is likely complex and these observations suggest that both the nanoscale organization of RyRs and the pattern of phosphorylated RyRs within clusters could be critical determinants of Ca2+ spark dynamics. To test this hypothesis, we used computational modeling to quantify the relationships between RyR cluster geometry, phosphorylation patterns, and sarcoplasmic reticulum (SR) Ca2+ release. We found that RyR cluster disruption results in a decrease in spark fidelity and longer sparks with a lower amplitude. Phosphorylation of some RyRs within the cluster can play a compensatory role, recovering healthy spark dynamics. Interestingly, our model predicts that such compensation is critically dependent on the phosphorylation pattern, as phosphorylation localized within the cluster center resulted in longer Ca2+ sparks and higher spark fidelity compared to a uniformly distributed phosphorylation pattern. Our results strongly suggest that both the phosphorylation pattern and nanoscale RyR reorganization are critical determinants of Ca2+ dynamics in HF.

Introduction HF is often characterized on the cellular level by a marked loss of t-tubules and the creation of "orphaned" RyRs. However, recent data from super-resolution imaging studies have revealed that the morphology of CRUs also changes in this condition, resulting in smaller and more numerous RyR clusters [4,[7][8][9] (Fig 1). This RyR cluster dispersion has been linked to slower Ca 2+ spark kinetics and a resulting desynchronization of the overall Ca 2+ transient [7,10]. However, RyR activity is also critically regulated by phosphorylation by protein kinase A (PKA), Ca 2+ calmodulin kinase type II (CaMKII), and various phosphatases [11] (Fig 1). It is well known that phosphorylation increases RyR open probability [11], and that RyR phosphorylation is augmented during HF [2,12,13]. Therefore, understanding RyR function and dysfunction during HF requires an integrated understanding of nanoscale RyR localization and phosphorylation status. Complicating the issue is the finding that RyR phosphorylation patterns may not be uniform. Indeed, Sheard et al. [8] reported while RyRs are uniformly phosphorylated across clusters in healthy cardiomyocytes, cells from failing hearts exhibited a higher density of PKA-phosphorylated RyRs at the center of RyR clusters. The functional implications of these changing phosphorylation patterns are unclear, and have not been addressed in previous computational studies which assumed that all RyRs in a cluster are equally phosphorylated [14][15][16].
In the present work, we investigated how the spatial organization of RyR clusters affects Ca 2+ dynamics, with a particular focus on changing patterns of phosphorylation. To this end, we have adapted an existing mathematical model of the CRU [7] to include distinct Ca 2+ sensitivities of individual RyRs. This is a theoretical computational study inspired by previous experiments with the goal of dissecting how phosphorylation of RyRs can impact Ca 2+ dynamics. We based the choice of our CRU geometries on a previous study from Kolstad et al. [7], where geometries were set based on CRU images from superresolustion microscopy data. For our current study, we adjusted those experimentally-based geometries in order to maintain consistency across simulations in terms of RyRs per geometry. The goal is not to replicate a specific set of experiments, but rather to understand how phosphorylation patterns may affect Ca 2+ kinetics, based on observations from previous experimental studies. We expect that our predictions will motivate the design of experiments that can decipher how these localized, nanoscale relationships contribute to impaired Ca 2+ homeostasis during HF.

Model development
The mathematical model from the work of Kolstad et al. [7] was adapted to differentiate between two subgroups of RyRs: phosphorylated and non phosphorylated RyRs. As described in [7], the model was extended from the model of Hake et al. [17] and includes a stochastic model of RyR opening developed by Cannell et al. [18] for both phosphorylated and unphosphorylated populations. The system of partial differential equations applied for the spatiotemporal evolution of [Ca 2+ ] in the cytosolic domain O c was: with c corresponding to the calcium concentration in the cytosolic domain, D c the diffusion constant of calcium, b i the concentration of the corresponding buffer and B i tot the corresponding total buffer concentration. In the cytosolic domain the buffers ATP, calmodulin, troponin and Fluo-4 were included and correspond respectively to numbers 1 to 4 in the functions listed in Eq (1). For the SR domain O s (both junctional and non-junctional SR), one calsequestrin buffer was included and its concentration will be denoted by b 5 with s corresponding to the calcium concentration in the SR domain. Both domains were coupled at the SR membrane with the following boundary conditions: The SERCA flux formulation was taken from the three-state SERCA model from Tran et al. [19]. The J RyR and J pRyR formulations were identical to the J RyR flux formulation from Kolstad et al. [7], where a two state stochastic model was used: Here, C describes the conductive state and O the non-conductive state. The k − and k + variables were defined as: The difference between J RyR and J pRyR is the value used for K + . The choice of the parameter K + for the phosphorylated and the non phosphorylated case is discussed in Section "K + parameter estimation". S1 and S2 Tables show the model and buffering parameters. The code and all relevant data can be found in the code repository.

Geometries
The model consists of a single CRU containing both cytosolic and SR domains (O c \ O s ). We assumed that the simulated CRU was a 3D cube with volume 1.008 μm × 1.008 μm × 1.008 μm. Cube dimensions were 12 nm × 12 nm × 12 nm (Δx = 12 nm) with 84 voxels per axis.
The RyR geometries were inspired by the super-resolution images of Kolstad et al. [7]. However, in this work we constrained the number of RyRs to 50 for all geometries, which is within the measured range in experiments [7]. Five different geometries were designed. The first two geometries (G1 and G2) contain a single cluster to simulate the healthy case (Fig 2A). The main difference between G1 and G2 is the shape of the cluster; while G1 is a more compact distribution, G2 is oblong. By comparing G1 and G2, we can determine whether altering cluster shape without adjusting the density has an impact on spark dynamics. The latter 3 geometries (G3, G4, and G5) differ in the number of CRU sub-clusters. Disrupted RyR clusters containing several sub-clusters have been observed in rats with HF [7]. Therefore G3, G4, and G5 contain 50 RyRs distributed in 2, 3, and 12 sub-clusters, respectively.
Each RyR has an area of 36 nm × 36 nm. The center-to-center distance between two neighbouring RyRs is 36 nm, as described in [7]. The jSR was designated as the area of a single receptor around each RyR ("padding"), in the same fashion as Kolstad et al. [7]. The SR surface area is given in S3 Table. The nSR was assumed to be a regular grid throughout the cytosol. The five different RyR geometries studied in the model. The red dots show the single RyRs and the tan coloured area represents the junctional SR membrane. All five geometries contain 50 RyRs. The first two geometries (G1 and G2) contain a single cluster to simulate the healthy case. The latter 3 geometries (G3, G4, and G5) differ in the number of CRU sub-clusters, to simulate the disrupted RyR clusters observed during HF. G3, G4 and G5 are organized into 2, 3 and 12 sub-clusters respectively. (B) The computational domain of the cytosol is presented. The space corresponds 1.008 μm x 1.008 μm x 1.008 μm. The grey bars represent the non-junctional SR. https://doi.org/10.1371/journal.pcbi.1010126.g002 Due to the diffusion model, both the Ca 2+ kinetics in the nSR and jSR are influenced by the net change in Ca 2+ concentration. The chosen structure matches the volume and surface area measured by Hake et al. [17]. The SERCA surface has a range of 4.54-4.85 μm 2 . The jSR and nSR surface areas agree with the surface areas measured by Hake et al. [17] in EM tomography dyad reconstructions. Fig 2B shows the location of the RyRs and jSR within the nSR grid for the G1 geometry.
The t-tubules where modeled creating a 1 voxel wide cleft space (12 nm) around the RyRs as described in [7].

Phosphorylation patterns
Since the goal of this study is to analyze the impact of spatial phosphorylation patterns on Ca 2+ spark signals, three different phosphorylation patterns were applied (Fig 3). Additionally, simulations with no phosphorylated RyRs and with blanket phosphorylation of all RyRs (to emulate previous computational models) were carried out. This resulted in 8 phosphorylation setups per geometry. An example of the 8 phosphorylation setups for geometry G1 is shown in Fig 3. Sheard et al. [8] reported that in HF, a spatial gradient of RyR hyperphosphorylation occurs with the highest phosphorylation levels occurring within the center of the nanodomain. In order to simulate this, phosphorylated RyRs were chosen using principal component analysis (PCA) [20]. The eigenvectors of the covariance matrix of the RyR positions were calculated to estimate the directions of maximal information of the CRU. With both orthogonal eigenvectors, an ellipse is generated. The ellipse dimensions were then decreased until the desired number of RyRs was phosphorylated. S1 Fig depicts an example of the PCA method applied to geometry G1 to calculate the inner phosphorylation. By applying this technique we select the desired number of RyRs in the center of the CRU when simulating HF conditions. We also chose to study the opposite pattern using the same PCA technique but phosphorylating RyRs near the outer boundary of the cluster. Although this specific phosphorylation pattern has not been observed in experimental studies, we include it in our simulations for completeness. These two PCA-based phosphorylation schemes will be referred to as inner and outer phosphorylation patterns from here onwards. For each pattern, two phosphorylation levels were simulated: 20% (10 RyRs) and 50% (25 RyRs) phosphorylation of the cluster. Finally, in order to simulate a uniform distribution of the phosphorylation pattern, as may be the case in the healthy condition (described in [8]), we simulated a uniformly distributed phosphorylation condition (again at 20% and 50% levels). In total, 8 different phosphorylation configurations per geometry were analyzed, all of which are visualized for the G1 geometry in Fig 3. An overview of the phosphorylation patterns used for G2, G3, G4 and G5 is provided in S2 and S3 Figs. Since the 8 different phosphorylation patterns were simulated for 5 different geometries, this results in 40 different spatial configurations.

K + parameter estimation
RyR phosphorylation increases the channel sensitivity to Ca 2+ [12]; we simulated this by lowering the K + value for a phosphorylated receptor (lower K + value for J pRyR than J RyR ). For the no phosphorylation condition, see red dots in Fig 3, the K + value was set to 55 μm. This value was chosen by optimizing K + values to achieve physiologically realistic spark property values [21]. A particular focus was set on obtaining non-phosphorylated sparks with time to peak (TTP) values of around 7 ms [21]. For phosphorylated RyRs, we differentiated between two specific cases. The first case is one where a single geometry contains both phosphorylated and non phosphorylated RyRs, see the outer, inner and uniform phoshorylation patterns in Fig 3. In this case, for phosphorylated RyRs (yellow dots in Fig 3) the K + value was set to 25 μm. This K + value was estimated based on obtaining sparks with spark properties within the same range as previous experimental results [6,21]. The selection of K + values for blanket phosphorylation was set such that the Ca 2+ dynamics from a blanket phosphorylation pattern would match the Ca 2+ kinetics obtained from using a uniform pattern. Computational studies generally assume all RyRs in a CRU to be identically sensitized [14][15][16]. Therefore, we wanted to compare this approach of a blanket phosphorylation assumption with assuming a uniform pattern, which is the pattern observed experimentally in healthy cardiomyocytes [8]. Since we assume two degrees of uniform phosphorylation (20% and 50%), we decided to use two different K + values for the blanket phosphorylation case, see Fig 3. K + was set to 45 to compare a blanket phosphorylation with a 20% uniform phosphorylation. Additionally, K + was set to 35 to compare a blanket phosphorylation with a 50% uniform phosphorylation. An overview of the chosen K + values is given in Table 1. The K + values were kept constant through all 5 different geometries, such that the only difference between the different geometries is the spatial localization of the RyRs.

Numerics
Since the model used assumes a stochastic model for RyR opening, 200 runs were carried out for each of the 40 configuration setups (5 geometries × 8 phosphorylation patterns). The simulations were initiated by randomly opening 1 RyR within the specific geometry, and the RyR fluxes were then calculated. Opening and closing of a single RyR was based on local Ca 2+ concentrations. The simulations were terminated when all RyRs were closed and remained closed for 1 ms. Therefore, we define spark duration as the duration of the simulations 1 ms after all RyRs are closed. Since the model is computationally intensive, simulations were generally terminated 1 ms after all RyRs, both phosphorylated and non phosphorylated, were closed, a point at which relevant metrics could be attained. After all RyRs are closed, the remainder of the spark is modeled as an exponential decay function which is independent of the RyR cluster spatial organization and the phosphorylation pattern. In S4 Fig, it can be observed that running the full simulation leads to successful termination of the Ca 2+ spark. The stochastic RyR model was calculated for a fixed time step of Δt = 1 ms. The model calculations were very stiff due to the small element volumes and the large fluxes. Therefore, the calculations of the fluxes in the model, described in the section "Model development", were solved analytically in the same fashion as [7]. We also used the following definitions to classify and quantify our results. We consider a spark to be successful if local increases in simulated Ca 2+ fluorescence intensity exceeded 30% (ΔF/F 0 � 0.3), and used the fraction of these successful sparks to estimate spark fidelity. The confidence interval of spark fidelity was calculated using the Agresti-Coull confidence interval [22]. Additional spark properties including the amplitude, time to peak (TTP), and spark simulation duration were used to compare and contrast different scenarios. For these three properties, the confidence interval was calculated using the standard error.

Circular CRUs increase spark fidelity and amplitude
First, simulations without RyR phosphorylation were carried out for all 5 geometries. S5 Fig shows that our simulation results match the results from Kolstad et al. [7], as spark fidelity and amplitude decrease with increasingly dispersed RyR configurations. We additionally compared two single-cluster geometries (G1 and G2), where G1 had a compact circular geometry and G2 a more oblong arrangement. Interestingly, a clear difference in spark fidelity (ie. proportion of sparks with ΔF/F 0 > 0.3) is observed between these two geometries (S5 Fig). This is illustrated by the spark time courses for G1 versus G2 (Fig 4A and 4B, respectively). For G1, 63 out of 200 sparks were successful (dark gray sparks), while for G2 only 29 successful sparks occurred. Fig 4C shows the probability of spark generation for both geometries. Out of 200 simulations, 31.8% ± 6.4% surpassed the spark detection threshold for G1. This number was significantly lower in the G2 geometry, which had a spark fidelity of 15.2% ± 4.9%. In addition to a decrease in fidelity, the oblong geometry in G2 also results in a significantly lower amplitude ( Fig 4D) and longer time to peak (TTP), see Fig 4E. A scatter plot of these measurements in individual sparks is shown in Fig 4F. We next investigated whether the lower amplitude, slower kinetics, and lower fidelity of Ca 2+ sparks generated by the oblong G2 vs compact G1 configuration were linked to a differing number of RyRs being activated. To this end, we plotted a histogram of the number of RyRs opened per simulation (Fig 4G). Note that the G1 geometry exhibits an "all-or-none" behavior, with bimodal clustering near either one or all 50 of the RyRs being activated in a single simulation. In comparison, G2 shows more varied behavior, and no simulations with full activation. This behavior is also observed in the activation map of a single simulation for both G1 and G2 (S6 Fig). For the G2 simulation activation map shown, the RyRs located at the edges of the geometry are not activated. Thus, our results show that the nanoscale organization of the unphosphorylated RyR clusters affects Ca 2+ dynamics, as compact, circular nanoclusters generate larger and faster Ca 2+ sparks than elongated nanoclusters.

Phosphorylation pattern of the RyR cluster determines spark properties
As shown in [8], different spatial patterns can be observed for RyR phosphorylation in the CRU. However, computational models generally assume all RyRs in a CRU to be identically sensitized due to a lack of spatial detail [14][15][16]. We next tested if there was a difference in Ca 2+ dynamics when all receptors are phosphorylated ("blanket" phosphorylation) as opposed to some percentage of uniformly phosphorylated receptors in the cluster (see Fig 3 for schematic). This allows us to determine whether the approach commonly used by modellers, assuming "blanket" phosphorylation, matches spark dynamics of uniformly distributed RyR phosphorylation observed in CRUs in healthy conditions [8]. To capture the effects of receptor phosphorylation and calcium sensitivity, we varied K + in the range [25 μm, 35 μm, 45 μm, and 55 μm]. Here 55 μm represents non phosphorylated RyR and therefore K þ RyR . On the other hand, 25 μm represents the situation of uniform phosphorylation (both 20% and 50%); 35 μm and 45 μm represent blanket phosphorylation scenarios. These three values represent K þ pRyR , see Table 1.
Plotting spark duration against spark amplitude (Fig 6B), we found that outer phosphorylation produces slightly larger sparks than uniform or inner phosphorylation (ΔF/F 0 = 0.560 ± 0.001 for 50% outer phosphorylation vs 0.547 ± 0.001 for 50% uniform phosphorylation vs 0.531 ± 0.002 for 50% inner phosphorylation, Fig 6C). We also measured TTP and spark duration for the same phosphorylation patterns (Fig 6D and 6E). Phosphorylation leads to increased spark duration in our model, across all patterns. An increased TTP was also seen across the different patterns, with the exception of the inner 20% configuration. For 20% phosphorylation, outer phosphorylation leads to the highest increase in TTP, followed by uniform, then inner ( Fig 6D). However, for 50% phosphorylation there is no significant difference between the TTP for the inner and uniform patterns (around 8 ms). This suggests that the relative effect on TTP decreases with increasing fraction of phosphorylated RyRs.
The sparks measured with an inner phosphorylation pattern have a longer duration compared to the uniform and outer case (15.2 ± 0.18 ms vs 13.4 ± 0.14 ms vs 13.0 ± 0.11 ms at 50% phosphorylation, and 13.0 ± 0.19 ms vs 12.4 ± 0.15 ms vs 12.7 ± 0.24 ms at 20% phosphorylation), despite the shorter TTP (see Fig 6F-6H). Importantly, phosphorylation patterns were found to have a large impact on the mean open time of individual RyRs within the cluster ( Fig  6I-6K). With a uniform phosphorylation pattern, RyRs located near the center of the cluster tended to be open longer than the channels located at the outer boundary (Fig 6J). In other words, the central channels anchored regenerative release. This effect is emphasized even more strongly with an inner phosphorylation pattern, since the inner channels are sensitized both by their spatial location and their phosphorylated state (Fig 6I). For the outer phosphorylation pattern, the effect is reversed, with channels across the cluster exhibiting close to uniform mean open times (Fig 6K). Later studies using super resolution imaging to visualize RyR CRUs have measured around 30 RyRs per cluster [6]. Therefore, we investigated whether the effects of phosphorylation pattern were also observed in clusters with 30 RyRs per cluster. S7 Fig shows that similar effects as in Fig 6 are observed when assuming a compact geometry containing 30 RyRs. For smaller geometries, the effect of the phosphorylation pattern in compact geometries is even more pronounced for fidelity, amplitude and TTP, while larger geometries lead to a higher impact on spark duration. Taken together, our simulations predict that the phosphorylation pattern within the CRU has a marked impact on Ca 2+ spark dynamics, as inner phosphorylation leads on average to higher fidelity, lower amplitudes, and longer spark durations than uniform or outer phosphorylation patterns.

Inner phosphorylation maximizes spark fidelity and increases spark amplitude in a dispersed RyR cluster organization
We next analyzed the effects of RyR phosphorylation in disrupted CRU geometries, which are reported in HF [7]. As noted above, when all the RyRs are unphosphorylated, we observed a step-wise decrease in spark fidelity and amplitude with increasing dispersal of the CRU geometries (S5 Fig). These results match the outcomes of previous work based on the same model [7]. As described by Sheard et al. [8], during HF, a shift of the PKA-mediated phosphorylation pattern is observed towards the center of the cluster. Therefore, we investigated whether the We first focus on the differences between G2 and G3, since the only difference between these geometries is that G2 contains one cluster and G3 two sub-clusters, while overall cluster shape is maintained (Fig 2). For both geometries, higher spark fidelity is observed when applying inner phosphorylation compared to a uniform or outer phosphorylation (Fig 7A). We also observe a significant increase in spark duration for both geometries with inner phosphorylation in comparison with uniform and outer phosphorylation patterns (14.2 ± 0.16 ms vs 12.9 ± 0.16 ms vs 12.76 ± 0.17 ms at major phosphorylation for G2, and 14.7 ± 0.24 ms vs 13.6 ± 0.22 ms vs 13.5 ± 0.24 ms at major phosphorylation for G3) (Fig 7B).
However, differing effects on spark amplitude were observed for compact geometry G2 and discontinuous geometry G3 (Fig 7C and 7D). Indeed, for G1 and G2, an increase in amplitude is observed when shifting phosphorylation from the inner to the outer pattern (Figs 6 and 7). G3 exhibits the opposite behavior, as the mean amplitude is highest when assuming inner phosphorylation on amplitude, followed by a uniform phosphorylation and an outer phosphorylation (0.532 ± 0.005 vs 0.500 ± 0.009 vs 0.412 ± 0.010 at major phosphorylation for G3). Interestingly, kernel density plots show that changing the phosphorylation pattern had rather complex effects on spark amplitude, as data from the G2 configuration showed a Gaussian relationship, while biphasic curves were observed for the G3 configuration. To understand the described effects of phosphorylation, histograms showing the number of opened RyRs vs number of simulations (Fig 7E and 7G) are illustrative. When plotting the no phosphorylation case for G2 and G3 in Fig 7E, we observe that in a successful spark for G3, only around 25 RyRs open, as activation was limited to one sub-cluster. For G2, however, we see that if a spark is generated the number of opened RyRs is much higher. It is notable that RyRs in the center of each subcluster tend to open more frequently and for longer durations than RyRs on the boundary (see Fig 7F) making it difficult for a spark to activate the neighbouring subcluster.
When applying a phosphorylation pattern, however, the activation map changed considerably (see Fig 7H). For the inner case, the number of open RyRs increased dramatically ( Fig  7G) and phosphorylated RyRs remained open for longer (Fig 7H), leading to more reliable activation of the neighbouring subcluster. This jumping of released Ca 2+ between subclusters was much less frequent with an outer phosphorylation pattern, with a lower number of RyRs being activated (Fig 7G).
To mechanistically understand the rather complicated dependence of spark kinetics on phosphorylation pattern, we studied the dynamics of individual simulations by tracking the opening and closing of the RyRs for single representative sparks (see S1 and S2 Videos). Based on these videos, we first confirm that with inner phosphorylation, the TTP is faster compared to the uniform phosphorylation pattern. However, due to the phosphorylation, the inner RyRs are more likely to reopen, which prolongs the total spark duration and creates a plateau of phosphorylation patterns in G1. (For no phosphorylation n = 63, for inner 20% n = 97, for inner 50% n = 125, for uniform 20% n = 94, for uniform 50% n = 121, for outer 20% n = 63, and for outer 50% n = 118). (D) TTP (mean ± standard error) for inner, uniform, and outer phosphorylation patterns in G1. (E) Spark duration (mean ± standard error) for inner, uniform, and outer phosphorylation patterns in G1.  Ca 2+ release, an effect that is particularly prominent at high phosphorylation levels. Additionally, the sensitized RyRs are able to sustain the regenerative release longer, resulting in higher fractional release from the SR and to allow the jumping of released Ca 2+ between subclusters.
The histogram also helps to explain the biphasic curves in the amplitude kernel density estimate plot (Fig 7D). The first peak, at lower amplitude, represents sparks generated when only one sub-cluster is activated. The second peak, at higher amplitudes, represents sparks generated when both sub-clusters are activated. For inner phosphorylation, the mode occurs at a higher amplitude of around 0.55, whereas outer phosphorylation more often reaches an amplitude of 0.35 more characteristic of single subcluster activation.
We next investigated the effect of phosphorylation on clusters with pronounced dispersion (Fig 8). The activation maps for G4 and G5 are shown in Fig 8A. We found that inner phosphorylation yielded higher spark fidelity pattern than for a uniform or outer phosphorylation pattern ( Fig 8B); effects that were more marked than observations made in condensed CRUs (compare with Fig 6). For 50% phosphorylation the fidelity is 38.7 ± 6.69% for an inner phosphorylation, whereas for uniform and outer phosphorylation, the fidelity drops to 30.4 ± 6.31% and 12.7 ± 4.57%, respectively. Indeed for G5, this effect is even more relevant since only the inner 50% phosphorylation and the uniform 50% phosphorylation patterns generate any sparks over the 0.3 threshold (Fig 8C). Spark properties calculated for G4 showed that, as for G3, the inner phosphorylation yields longer sparks with a higher amplitude ( Fig  8D-8F). We do not show the outcomes of the spark analysis for geometry G5 since the number of successful sparks is too low to generate meaningful statistics. These outcomes suggest that phosphorylation pattern may be of particular importance in severely disrupted CRU geometries, where reorganisation of phosphorylated RyRs toward cluster centers may provide a compensatory effect and allow sparks to persist.

Discussion
In this study, we have adapted a stochastic computational model to study the effect of RyR cluster organization and phosphorylation patterns on Ca 2+ spark dynamics. Using this model, we show that fundamental features of RyR nanoclusters in a CRU, including cluster geometry, phosphorylation pattern, and cluster integrity, interact in combination to regulate spark properties.
As shown in a previous computational study by Walker et al. [23], spark fidelity in randomly distributed non-square geometries is lower than fidelity in squared geometries. This matches with our observations regarding lower fidelity in less circular geometries, along with lower amplitudes and longer sparks. This can be understood by the fact that in a compact geometry, the RyRs are in closer proximity to other RyRs in the cluster, leading to tighter intra-cluster coupling of all RyRs within the cluster. Indeed, Ca 2+ leaks out of the sides of the CRU geometry into the cytosol, meaning that RyRs on the outer part of the CRU geometry can experience less Ca 2+ availability. In oblong clusters, such as the G2 geometry, these effects are amplified at the cluster ends where there are less neighbouring RyRs, leading to reduced estimate plot of the amplitudes for the three phosphorylation patterns in G3. (E) Histogram of opened RyRs per simulation for G2 and G3 with no phosphorylation. (F) Geometric visualization of mean open times for each RyR throughout all 200 simulations for G2 and G3 assuming no phosphorylation. (G) Histogram of opened RyRs per simulation for G3 with inner and outer phosphorylation. (H) Geometric visualization of mean open times for G3 with inner and outer phosphorylation. � =p<0.05, �� =p<0.01 and ��� =p<0.001. The p values were calculated with respect to the non phosphorylated case using a t-test, except for the fidelity p value which was calculated from Fisher's exact test (due to binary data).

PLOS COMPUTATIONAL BIOLOGY
coupling. To visualize this we have plotted the histogram of neighbours in S8 Fig. We observe that G1 RyRs generally have a higher number of neighbours than the RyRs in geometry G2.
We further observed that more dispersed CRU geometries result in decreased spark fidelity and amplitude, and increased spark duration, consistent with previous computational and experimental work [7]. These results and the differences between compact and disrupted geometries match the spark property changes measured between healthy and HF patients [24].
Our model adaptions have allowed the novel study of the effect of phosphorylation spatial patterns on Ca 2+ dynamics. Based on the redistribution of the spatial pattern observed using enhanced expansion microscopy for the RyR phosphorylation site pSer2808 on Wistar rats undergoing right ventricle HF [8], we studied the effect of spatial phosphorylation pattern redistribution on multiple cluster geometries. In our simulations, we observe that inner phosphorylation results in compact distribution of phosphorylated RyRs in the center of the cluster, whereas the outer and uniform patterns can be seen as disrupted phosphorylation subdomains. As seen in this work and in a previous study by Walker et al. [23], a compact cluster leads to higher fidelity, meaning increased likelihood for RyRs to open and release calcium. This increase in Ca 2+ makes neighbouring RyRs more likely to open. We find that the average amplitude of the spark is significantly higher when assuming an inner phosphorylation pattern than when assuming a uniform distribution or an outer phosphorylation pattern because there is more consistent activation of adjacent sub clusters. This can be observed in the representative S1 and S2 Videos, where inner phosphorylation leads to an activation of both subclusters, whereas the uniform phosphorylation only activates one of the subclusters. Cluster dispersion leads to low excitability of clusters, and imparied EC coupling fidelity, which is a common property of failing cardiomyocytes. Our simulations show that phosphorylation counteracts this by sensitizing channels and increasing excitability. Additionally, we have observed that the reorganization to an inner phosphorylation pattern can further strengthen this compensatory effect and is even more pronounced in disrupted geometries than in compact geometries.
The phosphorylation of the RyRs is a widely studied and controversial field, since several papers show contradictory results [25]. Marx et al. [12] report the dissociation of the FK506 binding proteins (FKBP12.6) from the RyRs after PKA phosphorylation. However, as reported by Bers [25], this finding doesn't match with experiments measured by other groups [26]. In a different study Marx et al. report the impact of FKBP12.6 on coupled gating between neighbouring RyRs [27]. Based on the results of Marx et al., Sobie et al. [28] used a computational model to study the impact of coupled gating between RyRs on calcium sparks. Their results justify the decrease in spark fidelity by a disruption of the coupling gating between RyRs. The works of Marx et al. and Sobie et al. were carried out before super resolution imaging techniques were applied to visualize single molecule localization (e.g. individual RyRs) [5]. Super resolution imaging has since proven the disruption of RyR clusters during HF [7]. In a recent study from Asghari et al. [29], they observe that phosphorylation of the RyRs also leads to disruption of cluster geometries. In their study, the disrupted geometries also show higher spark fidelity and longer full duration half maximum values. These results may seem to contradict the outcomes from Kolstad et al. [7], where disrupted geometries lead to lower fidelities and lower amplitudes. However, in the study from Asghari et al., healthy compact geometries are compared with phosphorylated disrupted geometries. In our study we also saw that disrupted geometries undergoing phosphorylation could lead to longer and more frequent sparks than compact geometries without phosphorylation (Fig 9). For example, G3 undergoing inner phosphorylation leads to higher fidelities and higher amplitude values than sparks generated from G1 unphosphorylated conditions. Thus, the results from this study are consistent with both the outcomes from Kolstad et al. [7], when assuming no phosphorylation, and with the results from Asghari et al. when comparing compact unphosphorylated geometries with dispersed phosphorylated geometries.
An additional controversial aspect is the impact of PKA versus CaMKII on Ca 2+ spark properties. Guo et al. [2] studied the effect of CaMKII on Ca 2+ sparks in mouse ventricular myocytes. They observe an increase in Ca 2+ spark frequency, spark duration, spatial spread, and amplitude due to CaMKII-mediated RyR phosphorylation. Furthermore, Li et al. [30] found that PKA increases Ca 2+ spark frequency, amplitude, duration, and width in mouse ventricular myocytes. However, they conclude that PKA mediated changes in spark dynamics may be attributable to phospholamban and its resulting effects on SERCA rather than RyR phosphorylation [30]. CaMKII and PKA therefore have different impacts on calcium sparks, although both phosphorylating proteins increase the open probability of the RyR [30]. Further computational studies investigating phosphorylation patterns for CaMKII-specific sites may be of interest to determine whether distinct phosphorylation patterns occur in CaMKII-mediated phosphorylation and to understand why CaMKII and PKA show different phosphorylation effects as reported by Guo et al. [2] and Li et al. For example, the larger size of CaMKII (56-58 kDa [31]) may hinder phosphorylation activity in the central part of the CRU in the narrow dyadic cleft, leading to an outer phosphorylation pattern.
In this study we reaffirm the importance of cluster geometry on Ca 2+ spark properties. Additionally, these results show the large impact that spatial phosphorylation pattern can have on Ca 2+ spark properties, leading to differences in spark fidelity, amplitude or duration depending on the pattern. It is necessary to perform further studies, both computational and experimental, to explore larger-scale impacts of RyR phosphorylation patterns on whole cell calcium dynamics. For example, a computational study including multiple nearby clusters may address whether specific phosphorylation patterns may allow for more reliable spread of calcium sparks into calcium waves through the dyadic region. Additionally, further simulations including an increase in the volume of the dyadic junction due to t-tubule loss in heart failure would be interesting. We believe that an increased volume should lead to a desensitization of the RyR channels to calcium and the effect observed during RyR disruption should be even more pronounced. Further imaging studies based on the spatial distribution of the phosphorylated RyRs as shown in [8] and on the relation between cluster geometry and phosphorylation during HF [7,29] should be considered. Future work using this model and other models may also allow us to define functional clusters using imaging data, as simulations can determine subcluster distances that allow for reliable spark production. Differences in cluster geometry and spatial phosphorylation patterns could potentially explain the conflicting results from different studies on the impact of phosphorylation proteins on Ca 2+ sparks [25], but also emphasize the need for future investigation of the role of PKA and CAMKII phosphorylation of RyRs during health and disease.