The Interplay of Rogue and Clustered Ryanodine Receptors Regulates Ca2+ Waves in Cardiac Myocytes

Ca2+ waves in cardiac myocytes can lead to arrhythmias owing to delayed after-depolarisations. Based on Ca2+ regulation from the junctional sarcoplasmic reticulum (JSR), a mathematical model was developed to investigate the interplay of clustered and rogue RyRs on Ca2+ waves. The model successfully reproduces Ca2+ waves in cardiac myocytes, which are in agreement with experimental results. A new wave propagation mode of “spark-diffusion-quark-spark” is put forward. It is found that rogue RyRs greatly increase the initiation of Ca2+ sparks, further contribute to the formation and propagation of Ca2+ waves when the free Ca2+ concentration in JSR lumen ([Ca2+]lumen) is higher than a threshold value of 0.7 mM. Computational results show an exponential increase in the velocity of Ca2+ waves with [Ca2+]lumen. In addition, more CRUs of rogue RyRs and Ca2+ release from rogue RyRs result in higher velocity and amplitude of Ca2+ waves. Distance between CRUs significantly affects the velocity of Ca2+ waves, but not the amplitude. This work could improve understanding the mechanism of Ca2+ waves in cardiac myocytes.


INTRODUCTION
Ca 2+ sparks due to the opening of clustered RyRs are the elementary Ca 2+ release events in normal cardiac myocytes (Cheng et al., 1993;Cheng and Lederer, 2008), which could occur in self-propagating succession along the length, and contribute to waves of elevated Ca 2+ concentration under some pathological conditions (López-López et al., 1995). Ca 2+ waves have been observed in a diversity of cells (Ridgway et al., 1977;Fabiato, 1983;Cornellbell and Finkbeiner, 1991) and studied experimentally and theoretically (Fabiato and Fabiato, 1972;Fabiato, 1985;Backx et al., 1989;Swietach et al., 2010). Generating Ca 2+ waves in myocytes is associated with RyRs gating and sarcoplasmic reticulum Ca 2+ overload (Petrovic et al., 2015;Williams et al., 2017). Quarky Ca 2+ release (QCR or Ca 2+ quark) with a small amplitude and a long duration arising from rogue RyRs is another significant Ca 2+ release mechanism (Wang et al., 2001;Cheng and Wang, 2002;Brochet et al., 2011;Shang et al., 2014). Hence, Ca 2+ waves are a natural consequence of regenerative Ca 2+ releases of both Ca 2+ sparks and quarks. There is, however, lack of studies to relate Ca 2+ waves to the interplay of Ca 2+ sparks and quarks from the junctional sarcoplasmic reticulum (JSR). Based on the Fickian diffusion of cytoplasmic Ca 2+ , a computational model was developed to show the effects of rogue RyRs on Ca 2+ waves under heart failure (Lu et al., 2010). Given the spark-width paradox from the Fickian diffusion models (Walker et al., 2014), the anomalous diffusion model can solve the problem and look more deeply into the mechanism of diffusion (Sato and Bers, 2011).
On the other hand, one of the challenges in developing models for Ca 2+ waves is the inconsistence between computational and experimental free Ca 2+ concentration in the cytoplasm ([Ca 2+ ] cyto ) (Izu et al., 2013). The computational results of [Ca 2+ ] cyto were ∼20 µM  under physiological conditions or even as high as ∼100 µM (Izu et al., 2001;Chen et al., 2014) under pathological conditions, which disagrees with the measured [Ca 2+ ] cyto of ∼1 µM (Williams et al., 1985;Takamatsu and Wier, 1990). Although a "wave front sensitization" model showed [Ca 2+ ] cyto of ∼1 µM (Keller et al., 2007), Sobie et al. indicated that elevated JSR Ca 2+ level is a critical factor to raise RyRs open probability (Sobie et al., 2017). Hence, JSR Ca 2+ regulation should be incorporated into the computational models of Ca 2+ waves to show the decrease of Ca 2+ flux through rogue and clustered RyRs as the JSR depletes (Sobie et al., 2004;Picht et al., 2011;Izu et al., 2013).
The objective of the study is to quantify the interplay of rogue and clustered RyRs on regulating Ca 2+ waves in cardiac myocytes. A two-dimensional (2D) model of Ca 2+ waves in the cytoplasm was proposed with considering the distribution of clustered and rogue RyRs on the JSR membrane. The anomalous subdiffusion of Ca 2+ in the cytoplasm and JSR Ca 2+ regulation were also included. The stochastic opening Ca 2+ release units (CRUs) of rogue and clustered RyRs was regulated by free Ca 2+ concentrations in both cytoplasm and JSR lumen. With these features, we showed the importance of rogue RyRs on the initiation and propagation of Ca 2+ waves.

Geometrical Model
Considering the quasi-isotropic diffusion of Ca 2+ in the cytoplasm (Izu et al., 2001), we adopted a 2D model to mimic Ca 2+ waves. Figure 1A shows the geometrical model of a cardiac myocyte, the x-and y-directions of which refer to the longitudinal axis and z-line, respectively. Baddeley et al. have experimentally observed the RyR distribution on the JSR membrane (Baddeley et al., 2010). Most RyR channels form regular arrays, defined as "clustered RyRs." Others are rogue RyRs uncoupled from the clustered RyRs. Clustered and rogue RyRs are randomly distributed. Figure 1B shows schematic representative of the CRU distribution on JSRs, which includes CRUs of clustered RyRs (∼22 RyR channels in a CRU) and CRUs of rogue RyRs (∼3 RyR channels in a CRU). CRUs of clustered RyRs (∼2 CRUs in a JSR) are surrounded by randomly distributed CRUs of rogue RyRs (∼8 CRUs in a JSR).

Governing Equations
Ca 2+ release events are simulated synchronously by a hybrid model, which consists of two parts: a model of Ca 2+ waves in the cytoplasm and a model of Ca 2+ blinks in JSRs. The reactiondiffusion system for Ca 2+ waves in the cytoplasm based on the anomalous subdiffusion model, including the distribution of clustered and rogue RyRs, is described as follows: where [Ca 2+ ] cyto is the free Ca 2+ concentration in the cytoplasm, t is time, x and y are the spatial coordinates, D x (=300 µm 2 s −1 ) and D y (=150 µm 2 s −1 ) denote the Ca 2+ diffusion coefficients for anisotropic diffusion. The anomalous subdiffusion order β is 2.25. J dye is the flux due to the Ca 2+ fluorescent indicator dye, Fluo-4-AM, in the cytoplasm. J buffer−cyto is the flux due to the endogenous stationary buffers. J pump is the pumping rate of SR Ca 2+ -ATPase. SR pumps will be started when [Ca 2+ ] cyto exceeds the resting Ca 2+ concentration level (0.1 µM). The detailed description is in the Appendix A in Supplementary Material. On the other hand, the balance equation for Ca 2+ blinks in each JSR is written as: where [Ca 2+ ] lumen is the free Ca 2+ concentration in the lumen of a JSR. J release−lumen denotes the Ca 2+ release flux caused by opening of clustered RyRs (J clustered ) and rogue RyRs (J rogue ). J buffer−lumen is the Ca 2+ flux due to the buffer, calsequestrin, in the JSR. J refill is the refilled Ca 2+ flux to the JSR. The detailed description is in the Appendix B in Supplementary Material. Various parameters of the dye and buffers in the cytoplasm and JSR lumen are listed in Table 1, similar to previous studies Kong et al., 2013).

Firing Probability of Rogue and Clustered RyRs
The firing probability per unit time for CRUs of rogue or clustered RyRs is determined by Ca 2+ concentrations in the cytoplasm and JSR (Györke and Gyorke, 1998;Qin et al., 2008Qin et al., , 2009, which can be expressed as: where P cyto and lumen refer to the firing probability per unit time of Ca 2+ release events controlled by [Ca 2+ ] cyto and [Ca 2+ ] lumen , respectively. The detailed description is in the Appendix C in Supplementary Material.

Numerical Solutions
The 2D computational domain of a cardiac cytoplasm (20 × 20 µm 2 ) was meshed with squares of 0.1 × 0.1 µm to simulate Ca 2+ release events from multiple JSRs. JSRs (i.e., yellow circles with radius of 0.3 µm in Figure 1A) are uniformly distributed in the computational domain with l x (2 µm) along x-axis and l y (0.8 µm) along y-axis. Moreover, CRUs are stochastically distributed at nodes within each JSR, which includes 8 CRUs  of rogue RyRs (N rogue = 8) and 2 CRUs of clustered RyRs (N clustered = 2). As shown in Figure 1C, Equations (1-3) were solved using a FORTRAN-developed program similar to a recent study (Chen et al., 2018). The shifted Grünwald formula of center difference (Tadjeran and Meerschaert, 2007) was used to discretize the fractional differential term in Equation (1) as: where g k = Ŵ(k−α) Ŵ(k+1) , Ŵ denotes the Gamma function. α= β−1 = 1.25, k is an integer with α <k< α+1, and h is the mesh size. Free Ca 2+ concentrations in the cytoplasm and JSR were calculated simultaneously. The variable time step algorithm was used. The zero-flux boundary condition was set to the 2D computational domain of a cardiac cytoplasm for the Monte Carlo simulations.

Interplay of Rogue and Clustered RyRs in Neighbor JSRs
We have recently studied the effects of rogue RyRs on single Ca 2+ sparks and quarks using the model in Equations (1-5), which was validated against the experimental measurements in rat cardiac myocytes (Chen et al., 2018). Here, we used the experimentallyvalidated numerical model to investigate whether a Ca 2+ spark could trigger rogue and clustered RyRs in neighbor JSRs or not. Snapshots of Ca 2+ release events in a computational domain of 20 × 20 µm 2 were taken at 10, 20, and 40 ms when a CRU of clustered RyRs was fired. Figures 2A,B show rogue RyRs at point (11.9, 9.6) and clustered RyRs at point (12.0, 9.6), respectively, activated by the Ca 2+ spark at point (10.0, 9.6). Moreover, the release of clustered RyRs could trigger other CRUs of clustered RyRs in neighbor JSRs with the help of rogue RyRs, as shown in Figure 2C. The results demonstrate that Ca 2+ sparks from the opening CRUs of clustered RyRs could activate CRUs of clustered RyRs and rogue RyRs in neighbor JSRs to increase [Ca 2+ ] cyto .

Initiation of Ca 2+ Waves
Ca 2+ waves could be triggered in a domain where [Ca 2+ ] cyto is higher than the resting Ca 2+ concentration level (Lu et al., 2010;Izu et al., 2013). A Ca 2+ spark with a large current or several neighbor Ca 2+ sparks could also trigger Ca 2+ waves (Chen et al., 2014). Here, several neighbor Ca 2+ sparks are used to activate Ca 2+ waves in the simulation. The lowest propagation velocity of Ca 2+ waves was detected in the range of 40-110 µm/s (Cheng et al., 1996). Since the disappearance of Ca 2+ waves is related to a progressive decline of the wave propagation velocity, the longitudinal velocity is assumed to be higher than 40 µm/s. Figure 3A shows the probability of inducing Ca 2+ waves as a  function of the number of Ca 2+ sparks initiating from a corner in a square of 20 × 20 µm 2 for 100 ms. The beginning level of [Ca 2+ ] lumen is 1.0 mM. The initial Ca 2+ sparks arise from one CRU of clustered RyRs. Since Ca 2+ quarks increase Ca 2+ concentration in the cytoplasm, they enhance the probability of inducing Ca 2+ waves in myocytes. Accordingly, Figure 3B shows the number of triggered Ca 2+ sparks in 100 ms. A single Ca 2+ spark could not form Ca 2+ waves while four neighbor Ca 2+ sparks with the help of rogue RyRs guarantee the formation of Ca 2+ waves.

Propagation of Ca 2+ Waves
Ca 2+ waves in Figure 4 were generated in a computational domain of 20 × 20 µm 2 and recorded at 10, 50, 100, and 150 ms from left to right when the beginning level of [Ca 2+ ] lumen is 1.0 mM. Four Ca 2+ sparks were initiated at points (18, 19.2), (18,  Figures 4A,B indicates much faster Ca 2+ waves and higher amplitude when the effects of rogue RyRs are included in the computational model. Moreover, we found the "sparkdiffusion-quark-spark" mode. Ca 2+ released from clustered RyRs diffuses to a neighbor JSR, rogue RyRs are firstly activated in a stochastic manner to form Ca 2+ quarks, and subsequently they make activation of clustered RyRs to produce a Ca 2+ spark. The CRUs on the next z-line repeat the process to release Ca 2+ in the cytoplasm. Figure 4C shows the computational results for line-scan images of Ca 2+ concentration with and without considering rogue RyRs at the line of y = 16.8. The results reveal that rogue RyRs accelerate the propagation of waves by triggering more Ca 2+ sparks. The longitudinal propagating velocity has mean ± SD value of 95.9 ± 8.0 µm/s, which agrees with the experimental records (typically 100 µm/s) (Takamatsu and Wier, 1990;Wier and Blatter, 1991). Furthermore, Ca 2+ concentration in the simulation is in the range of 0.1-3.8 µM. The computational predictions are close to the experimental measurements (0.5-1.2 µM; Williams et al., 1985), and less than previous simulations (Izu et al., 2001;Chen et al., 2013).

Importance of [Ca 2+ ] lumen -Dependent Regulation
In previous mathematical models of Ca 2+ waves, [Ca 2+ ] lumendependent regulation was not considered, and Ca 2+ fluxes from CRUs were related to the current and duration of Ca 2+ sparks and quarks only. We determined the effects of [Ca 2+ ] lumendependent regulation on properties of Ca 2+ waves when the beginning level of [Ca 2+ ] lumen was set to 0.3, 0.6, and 0.9 mM. Local variation in the waves is shown in Figure 5A. Ca 2+ waves are very sensitive to [Ca 2+ ] lumen . There is an exponential increase in the velocity of Ca 2+ waves with the increase of [Ca 2+ ] lumen , as shown in Figure 5B. Moreover, a threshold value of [Ca 2+ ] lumen (i.e., >0.7 mM) exists for generation of steady Ca 2+ waves. Figure 5C shows the amplitude of Ca 2+ waves linearly increase with [Ca 2+ ] lumen because of the large driving force ([Ca 2+ ] lumen -[Ca 2+ ] cyto ) of Ca 2+ sparks and quarks.

Effects of Parameters of Rogue RyRs
There are a large number of randomly-distributed rogue RyRs near clustered RyRs. Sensitivity analysis on Ca 2+ waves was performed with respect to the changes in CRU number of rogue RyRs in a JSR (N rogue ). Figure 6A shows computational results for line-scan images of [Ca 2+ ] cyto at the line of y = 16.8 when CRU number of rogue RyRs in a JSR are 2, 8, and 14. As shown in Figure 6B, the amplitude and velocity increase with the increased CRU number of rogue RyRs because of high Ca 2+ quarks and Ca 2+ sparks.
The Ca 2+ release per CRU of rogue RyRs is mainly determined by the current through rogue RyRs (I rogue ) and the duration of current flow (T rogue ). The current and duration of rogue RyRs are 0.15 pA and 20 ms in Figures 2-6. Table 2 presents properties of Ca 2+ waves for different values of (I rogue × T rogue ). There is a strong correlation between Ca 2+ release through rogue RyRs and wave properties. When the release time of rogue RyRs decreases by half to 10 ms with the current of 0.15 pA, the longitudinal velocity and amplitude of waves decrease. When the current increases from 0.15 to 0.3 pA with the duration of 20 ms, Ca 2+ waves have higher values of amplitude and longitudinal velocity. Hence, the release amount of rogue RyRs characterized by (I rogue × T rogue ) is a significant parameter to affect wave properties.
Lu et al. pointed out that rogue RyRs were scattered over the 2D plane randomly (Lu et al., 2010). However, fluorescent imaging showed that QCR events were detected almost at the same site as Ca 2+ sparks (Brochet et al., 2011). The stochastic gating of a cluster contains 10-100 RyRs in a JSR and the mean number of RyRs is ∼21.6 (Soeller et al., 2007;Baddeley et al., 2009). Figure 7 shows that the distance between CRUs in the range of 0.05-0.2 µm has slight effects on the amplitude of Ca 2+ waves. However, the longitudinal velocities are 144.3 ± 13.1 and 63.6 ± 7.8 µm/s, respectively with respect to the distance of 0.05 and 0.2 µm. Therefore, the distance between CRUs of rogue RyRs and clustered RyRs should be taken into consideration in the propagation of Ca 2+ waves.

Effects on Activities From Subcellular to Cellar Levels
When a myocyte is under paced, several spontaneous Ca 2+ sparks or quarks can occur in the cytoplasm. It is, however, difficult to induce Ca 2+ waves that require abundant currents to trigger cardiac action potentials (Izu et al., 2013). In some  pathological conditions (e.g., arrhythmias), Ca 2+ waves occur in cardiac myocytes and affect the heart's normal function (Lakatta and Guarnieri, 1993). Ten Tusscher and Panfilov developed a ventricular cell model including subspace calcium dynamics that controls L-type calcium current and calcium-induced calcium release (CICR) (ten Tusscher and Panfilov, 2006). The model was used to study the effects of I Na recovery dynamics in combination with action potential duration (ADP) restitution on alternans and spiral breakup. On the other hand, the present results suggest that Ca 2+ release from JSRs is prone to be initiated to cause Ca 2+ waves with the help of rogue RyRs when [Ca 2+ ] lumen is large enough. Moreover, when the CRU number of rogue RyRs in a JSR or Ca 2+ release per CRU of rogue RyRs is large enough, or the distance between CRUs is small, the membrane potential would be elevated significantly to trigger action potential in single myocytes.

Potential Implications
This study shows some implications to heart diseases relevant to Ca 2+ waves in cardiac myocytes. According to our simulations, JSR Ca 2+ overload could increase RyR opening probability and generate Ca 2+ waves in heart . Mutation in RyRs could trigger ventricular tachycardia and sudden cardiac death (Lehnart et al., 2006). Moreover, the amplitude and velocity of Ca 2+ waves are significantly affected by the parameters of rogue RyRs, which may contribute to the formation of fibrillation (Macquaide et al., 2015) and arrhythmias (Ter Keurs and Boyden, 2007). A reduction of CRU number, the averaged current and releasing time of rogue RyRs resulted in an inhibition of Ca 2+ waves or dyssynchronous Ca 2+ transients in myocytes of congestive heart failure (Louch et al., 2013). Therefore, inhibition of Ca 2+ quarks through rogue RyRs may be a promising therapeutic target to prevent fibrillation in congestive heart failure.

Critique of the Study
The present study showed Ca 2+ waves relevant to the interplay of rogue and clustered RyRs. It addressed the importance of rogue RyRs to increase the initiation of Ca 2+ sparks, the incidence and propagation of Ca 2+ waves when [Ca 2+ ] lumen is large. The amplitude and velocity of Ca 2+ waves are in agreement with the experimental measurements. However, the model should be improved such that more parameters are added to investigate the mechanisms of Ca 2+ waves. This study used the fixed release time of rogue and clustered RyRs. It should be a random variable depending on JSR regulation. The detailed structure of clustered RyRs should be taken into consideration because it can influence the frequency of Ca 2+ sparks (Walker et al., 2014). Moreover, the present study comes from the assumption of a 2D model in healthy myocytes. A 3D model should be developed to investigate Ca 2+ waves in future studies.

CONCLUSION
We developed a mathematical model to investigate the interplay of rogue and clustered RyRs on regulating Ca 2+ waves in the cytoplasm. The computational results on Ca 2+ waves agree with experimental measurements in cardiac myocytes. It shows that four neighbor Ca 2+ sparks at the corner of a cardiac myocyte could induce Ca 2+ waves. Ca 2+ quarks increase the probabilities of triggering Ca 2+ sparks and speed up the propagation of Ca 2+ waves at high [Ca 2+ ] lumen . A new wave propagation mode of "spark-diffusion-quark-spark" is put forward. Particularly, Ca 2+ waves could occur only when [Ca 2+ ] lumen is higher than a threshold value of 0.7 mM. More rogue RyRs in a JSR result in more opening CRUs of rogue and clustered RyRs. Besides, Ca 2+ release from CRUs of rogue RyRs is a strong factor of wave properties. The velocity of Ca 2+ waves is affected significantly by the distance between CRUs. This study helps to understand basic mechanisms of Ca 2+ waves in cardiac myocytes.

AUTHOR CONTRIBUTIONS
XC, YH, and WT designed and performed the numerical calculation; YF analyzed and interpreted data; XC, YH, and WT wrote the manuscript.

ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China (11732001 and 11328201) and the Leading Talents of Guangdong Province Program. We would like to thank Xi Chen for valuable discussions about simulation and revision of the manuscript.