Disorder in Ca²⁺ release unit locations confers robustness but cuts flexibility of heart pacemaking

Each excitation–contraction coupling cycle in the heart begins with the generation of rhythmic excitation in the sinoatrial node (SAN). To satisfy a given blood supply demand, cardiac muscle performance, defined by its state of contractile apparatus, Ca cycling proteins, and cell excitability must be in balance with the rate and rhythm of the excitation impulses generated by the SAN. Thus, the cardiac pacemaker function is a vital part of the excitation–contraction coupling as it “sets the stage” for timely interactions for all further downstream mechanisms.

In turn, the generation of normal rhythmic cardiac impulses is executed via timely interactions within and among SAN pacemaker cells that involve coupled signaling of both cell membrane ion channels and Ca cycling, dubbed a coupled-clock system (Maltsev and Lakatta, 2009; Lakatta et al., 2010). A key element of the system is the sarcoplasmic reticulum (SR) that rhythmically generates diastolic local Ca releases (LCRs) via Ca release channels, ryanodine receptors (RyRs). The LCRs contribute to SAN cell pacemaker function via activation of Na/Ca exchanger (NCX) current (I_NCX) that accelerates the diastolic depolarization (Huser et al., 2000; Bogdanov et al., 2001; Lakatta et al., 2010). Initially, the role of Ca release in cardiac pacemaker function has been numerically studied in the so-called “common pool” models (Kurata et al., 2002; Himeno et al., 2008; Maltsev and Lakatta, 2009; Imtiaz et al., 2010; Severi et al., 2012), in which spatial Ca dynamics was neglected and cell Ca release was approximated by a single variable (i.e., providing nonspatial, whole-cell description).

While such traditional simplified approach has yielded substantial progress in our understanding of Ca release role in pacemaker function (including modern theories of coupled-clock function [Maltsev and Lakatta, 2009] and ignition [Lyashkov et al., 2018]), it has fundamental limitations, and further progress requires a new type of modeling that would take into account spatial interactions (Maltsev et al., 2014). A general issue is that we still do not have a clear quantitative understanding of how the Ca clock and the coupled clock system emerge from the scale of molecules (such as RyRs) towards the whole-cell function (Weiss and Qu, 2020). Thus, how spatially distributed and whole-cell signaling events that coordinate cell function emerge from stochastic openings of individual release channels remains an unresolved fundamental problem despite extensive investigation.

A concrete issue with common pool models is that they operate with Ca signals in sub-micromolar range, that is, at least two orders of magnitude lower than in the local vicinity of the RyR release clusters, i.e., inside Ca sparks and diastolic LCRs (up to tens and even hundreds of micromolars [Stern et al., 2013; Stern et al., 2014]). The molecules involved in pacemaker function (ion channels, exchangers, pumps, and Ca-sensitive enzymes), however, operate via the high local Ca concentrations rather than a whole-cell average. Common pool models thus are analogous to a mean-field theory approach in physics, which usually does recover faithfully qualitative information about system behavior, but fail to zone in on the critical parameter values.

While the membrane clock operates as a limit cycle oscillator (Kurata et al., 2012), the Ca clock seems to operate by completely different mechanisms based on phase transitions (Maltsev et al., 2011) or criticality (Nivala et al., 2012; Weiss and Qu, 2020). RyRs are organized and operate in clusters of 10–150 channels (Greiser et al., 2020), known as Ca release units (CRUs) in all cardiac cells, including SAN cells (Maltsev et al., 2011; Stern et al., 2014). Ca release in ventricular myocytes occurs mainly via Ca sparks (defined as a single CRU Ca release) tightly controlled by L-type Ca channel (LCC) openings (local control theory; Stern et al., 1997), whereas LCRs in SAN cells consist of multiple sparks that emerge spontaneously via self-organization by means of positive feedback provided by Ca-induced Ca release (CICR). The synchronized CRU activation leads to oscillatory, phase-like transitions in SAN cells (Maltsev et al., 2011) that generate a net diastolic LCR signal that ultimately interacts with NCX and LCCs (Maltsev et al., 2013; Lyashkov et al., 2018) to ignite pacemaker action potentials (APs), which comprises the ignition theory of APs (Lyashkov et al., 2018). The probability that a Ca spark can “jump” in this way to activate its neighbor and form a propagating multispark LCR instead of remaining an isolated event depends on various parameters, including the amplitude of the spark (sometimes called I_spark; Zhou et al., 2009; Maltsev et al., 2011) and its nearest neighbor distance (Stern et al., 2014). The difficult and fascinating phenomenon that occurs in this system is that the area affected by propagation of an LCR (similar to connected component in percolation theory) is discontinuous in the probability of the Ca spark jump. This discontinuity is a phase-like transition first demonstrated in Maltsev et al. (2011), and the exact parameter values at which it occurs is known as criticality. Thus, to understand the intrinsic mechanisms of SAN cell operation, it is important to quantitatively predict critical parameter values at which the system changes its operational paradigm from sparks to LCRs.

Initially, the LCRs were numerically studied with CRUs located in a perfect rectangular lattice (Maltsev et al., 2011; Maltsev et al., 2013). These studies showed that the ability of an LCR to propagate regulates the size and the impact of the LCRs on the diastolic depolarization. However, RyR immunofluorescence exhibited notable disorder in the CRU distribution (Stern et al., 2014; Maltsev et al., 2016) and a more recent numerical model (Stern et al., 2014) was developed with two fixed sizes of CRUs. The larger CRUs were located in a rectangular lattice that lacked release propagation until smaller clusters were introduced in the CRU network providing bridges for release propagation among the larger CRUs. While this CRU geometry is closer to realistic distribution and provided functionality to the SAN cell model, the total number of CRUs was not controlled, and thus whether or not disorder in the CRU distribution per se is capable of impacting SAN cell function has not been examined and represents the specific aim of the present study.

Here, we approach the problem via numerical model simulations of SAN cell function with the same number of identical CRUs, but different spatial CRU distributions. Specifically, we test how various degrees of disorder (or noise) in CRU positions would influence SAN cell operation. We found that the least robust SAN cell function is achieved when CRU positions are distributed in the perfect rectangular lattice. As disorder in the CRU positions increases, nearest neighbor distances decrease, creating “shortcuts” for Ca release propagation resulting in AP firing rate increase and robust pacemaker function. The most robust function was achieved when CRUs were distributed independently uniformly randomly.

However, the robust function comes at the cost of the chronotropic reserve of the CRU-based mechanism: disorder-facilitated broad propagation of Ca release in the basal state recruits a major fraction of CRUs, leaving only a minor fraction of CRUs unrecruited. And vice versa, in the case of square lattice CRU distribution with longer nearest-neighbor distances, a larger fraction of CRUs remains ready to fire and not recruited in the basal state (creating the “reserve”). This reserve becomes available to accelerate AP firing during β-adrenergic receptor (βAR) stimulation. Thus, order/disorder in spatial CRU distribution provides a novel subcellular mechanism of cardiac pacemaker regulation that a SAN cell may utilize to achieve a perfect balance between robustness and flexibility commensurate with its specific location and function within SAN tissue.

In addition to self-organization of CRU firing via propagating CICR, CRUs can be recruited to fire via direct activation of LCC. Studies in transgenic mice have demonstrated specific importance of the Ca_v1.3 LCC isoform for cardiac pacemaker function (Mesirca et al., 2015). While both Ca_v1.2 and Ca_v1.3 isoforms are colocalized with RyR clusters (i.e., CRUs; Christel et al., 2012), Ca_v1.3 has a lower (more hyperpolarized) voltage activation threshold, thereby providing a larger contribution to diastolic depolarization by generating an inward current and by regulating RyR-dependent Ca release during SAN pacemaker activity (Torrente et al., 2016). Thus, the present study also tested possible importance of Ca_v1.3 for rabbit SAN function with respect to its interplay with CRU distribution.

Common pool models, such as those of Kurata et al. (2002) or Maltsev and Lakatta (2009), are on the cellular scale, lumping the contributions of all CRUs into one variable describing Ca release and one pool of junctional SR (JSR) with fixed total volume. On the other hand, a 3-D model of SAN cell developed by Stern et al. (2014) featuring interactions of individual RyRs provides too many fine details of intra-CRU Ca dynamics at a much higher computational cost and which are not so important for understanding Ca dynamics at the level of CRU to CRU.

Thus, an adequate model for our study should deal with Ca release at the submicron level (as CRUs detected by confocal microscopy). We have previously developed such a model (submicron, CRU-resolved model) featuring a 2-D array of stochastic, diffusively coupled CRUs located under cell membrane (Maltsev et al., 2011), and then added the full electrophysiological membrane clock (Maltsev et al., 2013). While approximation of Ca release at the submicron level suits our study, here we performed a further major model update, including a model of Ca clocking and an upgrade to the 3-D model (Fig. 1). The full details of the updated model are given in the Appendix, but here is a summary of important changes.

(1) A major weakness of our old model was that it approximated the CRU as having a built-in “refractory period” when the CRU could not activate. Thus, it was not a mechanistic “coupled-clock” model because the refractory period of Ca release (being the essence of the clock) was in fact an independent model parameter obtained as a direct read-off value from experimental data (Maltsev et al., 2011). After the refractory period, the CRU was placed into a “ready” state, meaning that it is allowed to open with a given dyadic cleft-Ca–dependent probability per unit time. Here, we modified the model to predict the timing of spark activation based on JSR Ca concentration using the present knowledge in this research area. A typical spark generated by our new model is shown in Fig. S1. CRUs are placed in the ready state when the JSR Ca reaches a certain threshold level of Ca as it refills during diastole via SERCA pumping and intra-SR diffusion. The threshold existence and precise value are suggested by experimental studies (Vinogradova et al., 2010), numerical model simulations via common pool models (Maltsev and Lakatta, 2009; Imtiaz et al., 2010), RyR-based models (Stern et al., 2014), and by a recent theoretical study based on Ising formalism (Veron et al., 2021).

(2) In the previous model, CRU Ca release flux was fixed to certain I_spark value (another independent model variable). In the new model, it is proportional to difference in [Ca] inside and outside the JSR.

(3) In our prior model, the stochastic CRU firing generated a spark of a fixed duration that was another independent model variable tuned to experimental data. In our new model, Ca release of a CRU is terminated when I_spark reaches a small value that is comparable with the release amplitude of a single RyR at a given SR Ca level, in accord with the spark death mechanism via disruption of inter-RyR CICR. This is in line with modern concepts of spark termination such as induction decay (Laver et al., 2013), spark “death” (Stern et al., 2013), and Ising formalism (Maltsev et al., 2017a). The basic idea of all three conceptions is that a Ca spark is generated and kept alive by positive feedback of inter-RyR CICR, until the SR becomes sufficiently depleted for RyR currents to wane and for inter-RyR CICR to be disrupted, culminating in Ca spark termination.

(4) The new model features a more realistic Ca distribution within the cell (Fig. 1). Our previous model was a 2-D model: CRUs resided in a square lattice under the cell surface with common pools of SR and cytosol; local (spatially resolved) Ca concentrations were predicted by the model only under cell membrane in 2-D. The new model includes three layers of intracellular voxels. Thus, it becomes a 3-D model, making cytoplasmic Ca dynamics more realistic and precise. Its structure is similar to that of the Stern et al. (2014) model, but the number of voxel layers from cell surface to cell center is limited to three (for more computational efficiency): fine submembrane voxels, intermediate “ring” voxels, and a larger cell core voxel. Importantly, the voxel layers are introduced only for computational efficiency and do not have any physical barriers or gradients. The submembrane voxels (i.e., submembrane space) are not artificially separated from the rest of the cell. All intracellular voxels, including those of submembrane space, have the same characteristics of Ca diffusion and buffering.

(5) We introduced JSR connected to free SR (FSR) via a diffusional resistance that determines the kinetics of JSR refilling with Ca. This is a key part of the Ca clock mechanism (Maltsev and Lakatta, 2009; Imtiaz et al., 2010; Vinogradova et al., 2010).

(6) Ca in JSR is buffered by calsequestrin.

(7) The CRUs in the prior model positioned as a square lattice. Here, we added the capability to introduce various degrees of disorder around the ideal square lattice positions of the CRUs so that the resulting distribution of nearest-neighbor distances would be spread and thus more closely reproduce the CRU network geometry reported in confocal microscopy studies (Stern et al., 2014).

(8) I_CaL formulation was improved by including a contribution of Ca_v1.3 isoform to I_CaL in addition to the cardiac-specific isoform Ca_v1.2. Both Ca_v1.2 and Ca_v1.3 isoforms are localized in CRUs (Christel et al., 2012), but have different voltage activation thresholds (see Appendix). Our simulations were initially performed with I_CaL comprised of cardiac-specific isoform Ca_v1.2. Then, we investigated effects of inclusion of Ca_v1.3 isoform into I_CaL (see the last two sections of the Results).

Fig. S1 shows typical sparks generated by our new model. Fig. S2 shows that SAN cell model with uniformly random distribution of CRUs features more robust spontaneous AP firing (right) versus that with a square lattice distribution of CRUs (left). Fig. S3 shows that the effect of βAR stimulation wanes as disorder in CRU positions increases. Fig. S4 shows that Ca_v1.3 accelerates AP firing rate by interacting with CRUs and I_NCX in the model with square lattice distribution of CRUs. Fig. S5 shows that Ca_v1.3 accelerates AP firing rate by interacting with CRUs and I_NCX in the model with uniformly random distribution of CRUs. Fig. S6 shows that Ca_v1.3 increases the robustness of pacemaker function. Fig. S7 shows the importance of LCC coupling to CRUs. Video 1 shows Ca dynamics in a SAN cell model with square lattice distribution of CRUs. Video 2 shows Ca dynamics in a SAN cell model with perturbed square lattice distribution of CRUs with SD = 0.25 μm. Video 3 shows Ca dynamics in a SAN cell model with perturbed square lattice distribution of CRUs with SD = 0.5 μm. Video 4 shows Ca dynamics in a SAN cell model with perturbed square lattice distribution of CRUs with SD = 0.75 μm. Video 5 shows Ca dynamics in a SAN cell model with uniformly random distribution of CRUs. Video 6 shows stochastic individual sparks generated by a dormant cell model with square lattice distribution of CRUs. Video 7 shows subthreshold oscillatory LCR signals generated by a dormant cell model with uniformly random distribution of CRUs. Data S1 provides the computer code used in this paper.

First, we tested how SAN cell function changes if the same number of identical CRUs is distributed differently under cell membrane. We generated and tested cell models having five different types of spatial distributions of CRUs with gradually increasing disorder in the CRU positions. Left panels in Fig. 2 show the cell cylinder surfaces “unwrapping” to squares with different simulated CRU distributions: (1) CRUs placed exactly at the nodes of a square lattice of 1.44 µm size. (2) CRUs slightly deviating from the “square lattice” positions, following Gaussian distribution with SD = 0.25 µm. (3) CRUs moderately deviating from square lattice positions, following Gaussian distribution with SD = 0.5 µm. (4) CRUs strongly deviating from the square lattice positions, following Gaussian distribution with SD = 0.75 µm. (5) Uniformly independently random CRU positions excluding overlap.

As disorder in CRU position increased, the SD of CRU nearest neighbor distance distribution increased, but the mean values decreased from 1.44 to 0.777 µm, respectively (Fig. 2, right). For a comparison, confocal microscopy measurements performed previously in six rabbit SAN cells reported average nearest distances between CRUs in each cell ranging from 0.71 to 0.89 μm (Stern et al., 2014), i.e., close to our simulated values in uniformly random and strongly perturbed CRU settings. In all these CRU settings, our numerical model simulations (with 100% Ca_v1.2 in I_CaL) generated rhythmic spontaneous AP firing. However, the rate of the spontaneous firing substantially increased as disorder degree increased and average AP cycle length (CL) shortened (Fig. 3; for specific values, see Table 1, Basal state row). The effect of randomness was substantial: in the extreme case of uniformly random distribution, the diastolic depolarization duration was shortened approximately in half (Fig. 4 A, red versus black traces), and the CL reduced by >30%, accordingly. This effect is comparable or even exceeds the effect of βAR stimulation on CL reported in rabbit SAN cells (Vinogradova et al., 2002). The average CL exhibits a linear dependence versus average nearest-neighbor distance (R² = 0.96; Fig. 3 A, inset), indicating importance of CRU interactions with their nearest neighbors (i.e., CICR) in this phenomenon.

Because our hypothesis was that spatial CRU distribution affects CICR among the CRUs, we examined the dynamics of CRU recruitment in our five RyR settings by monitoring N_CRU, i.e., the number of CRU firing at each moment of time during diastole. To compare dynamics of CRU recruitment during diastolic depolarization, we overlapped the traces of representative cycles for each CRU setting and synchronized them at the maximum diastolic potential (MDP; Fig. 4, A and B). We found that for all CRU settings, a major fraction of CRU pool became recruited to fire during diastolic depolarization, i.e., before the AP upstroke. However, uniform random CRU distribution resulted in a more synchronized and much faster recruitment of CRU firing compared to square lattice setting, with increasing degrees of disorder yielding an increase of recruitment rate.

As a result of faster CRU recruitment, the I_NCX and I_CaL were activated earlier in the cycle, accelerating diastolic depolarization and thereby shortening the CL (Fig. 4, C and D). In all scenarios, I_CaL was activated concurrently with I_NCX or after a short delay (Fig. 5), indicating that the system operates via the coupled-clock paradigm (Maltsev and Lakatta, 2009) and in accord with ignition theory (Lyashkov et al., 2018). The short delay increased with the increasing disorder in CRU positions, indicating that enhanced CRU recruitment activates I_NCX at low voltages, below I_CaL activation (less than −50 mV). To illustrate the diastolic CRU recruitment, we made representative screenshots of CRU firing for each CRU setting overlapped with instant local Ca distribution under the cell membrane at the membrane potential (V_m) of −45 mV that is close to the AP ignition onset and I_CaL activation threshold (Fig. 6). Mainly individual sparks were observed in the square lattice setting. In contrast, the uniform random setting showed propagating LCRs. The intermediate setting with various degrees of disorder exhibited intermediate recruitment. Different CRU firing patterns can also be seen in the respective videos of local Ca dynamics in subspace (Videos 1, 2, 3, 4, and 5).

SAN cells are very heterogeneous in their biophysical properties and their key ion current densities vary substantially, e.g., up to an order of magnitude for I_CaL (Honjo et al., 1996; Monfredi et al., 2018; Fig. 7 A, green circles). On the other hand, our simulations showed that the rate of spontaneous AP firing can be substantially reduced because of differences in spatial CRU distribution (Fig. 3). We then hypothesized that it is possible to find a range for realistic I_CaL densities at which spontaneous AP firing is impossible for the square lattice CRU distribution, but possible for uniformly random model. To match realistic I_CaL density values (in pA/pF) and our model parameter g_CaL (in nS/pF) determining maximum I_CaL conductance, we performed voltage-clamp simulations with the voltage step protocol similar to that in experimental studies (Honjo et al., 1996; Monfredi et al., 2018). To instantly attain steady-state activation and inactivation gating, we set respective gating variables fL and dL in the model to their steady-state values calculated at the holding potential of −45 mV. Thus, our basal state conductance 0.464 nS/pA in the model generated the maximum peak I_CaL current of 11.17 pA/pF shown by upper edge of the magenta band in Fig. 7 A that crosses the range of the experimentally measured I_CaL densities approximately in the middle. Please note that the graph in Fig. 7 A has a dual y-axis scale for I_CaL and g_CaL, respectively, bridging experimental and simulated data.

Then, we performed a wide-range sensitivity analysis for g_CaL from its basal value of 0.464 nS/pF (100%) down to 0.2552 nS/pF (55%) that still remained within the range of realistic I_CaL values (magenta band in Fig. 7 A). The analysis revealed a wide margin of g_CaL values (shown by red shade in Fig. 7 B) where the AP firing is impossible with square lattice distributions of CRUs (red circles), but possible with uniformly random CRU positions. All simulated traces of 16.5 s duration are shown in Fig. S2, and an example of the result is shown in Fig. 7 C for g_CaL of 0.3248 nS/pF.

This analysis revealed another important aspect of SAN function with different CRU distributions: the model with the square lattice CRU distribution was capable of generating AP firing with relatively long CL of 536 ms on average (before it failed), whereas the uniformly random CRU model could reach only 438 ms. Thus, SAN cells can (in theory) harness their CRU distribution to safely reach low rates.

A newly discovered paradigm of pacemaker cell function is the ability of some SAN cells to reversibly switch to dormant (nonfiring) state. It was found in isolated cells (Kim et al., 2018; Tsutsui et al., 2018; Tsutsui et al., 2021) and in SAN tissue (Bychkov et al., 2020; Fenske et al., 2020). Furthermore, nonfiring cells in SAN tissue can generate subthreshold signals that may be important signaling events for generation of synchronized cardiac impulses (Bychkov et al., 2020). Thus, we hypothesized that CRU distribution is important not only for AP firing, but also for the generation of subthreshold signals in dormant cells. One important factor in cell dormancy is a lower density of I_CaL (Tsutsui et al., 2021). Using the results of our sensitivity analysis, we chose a low g_CaL of 0.2552 nS/pF, at which AP firing is impossible in both uniformly random and square lattice models and compared subthreshold V_m signals in the nonfiring cells in each case (Fig. 8 A). We found that subthreshold V_m oscillations in case of uniformly random distribution are much more powerful and occurred at a faster frequency than those in the square lattice model (Fig. 8 A, inset, Videos 6 and 7) as evident from the overlapped power spectra of both V_m signals computed for the time period after AP failure (Fig. 8 B).

Normal pacemaker cell function is not only robust, but also flexible, i.e., able to increase the AP firing rate under stress (fight-or-flight reflex) or decrease it at rest. We tested if CRU distribution plays a role in the flexibility of SAN cell function by examining responses of our model to βAR stimulation. Our simulations showed that all five CRU settings from square lattice to uniformly random exhibit fight-or-flight reflex. The CL notably decreased in the presence of βAR stimulation (red bars versus blue bars in Fig. 9 A). All intervalograms are shown in Fig. S3, and examples of βAR stimulation effect in extreme cases of CRU distributions are shown in Fig. 9 B. The presence of disorder substantially decreased the stimulation effect; that is, the CL shortening decreased as disorder in CRU positions increased (both absolute and relative changes are given in Table 1). Furthermore, the absolute decrease in the CL was linearly related (R² = 0.9931) to the basal CL before βAR stimulation (Fig. 9 C) that is in accord with a recent experimental report (Kim et al., 2021).

To get further insights into the effect of spatial disorder on βAR stimulation, we simulated and compared the dynamics of N_CRU and I_NCX during a representative cycle for square lattice and uniform models (Fig. 10). We found that in both cases the βAR stimulation effect, linked to earlier and stronger recruitment of CRUs, translated to respective earlier and stronger activation of I_NCX. However, this effect (i.e., the shift to earlier CRU recruitment) was stronger in square lattice setting, indicating the presence of a larger reserve CRU pool ready to be activated and synchronized via βAR stimulation.

Thus far, we reported our results about functional importance of CRU distribution with I_CaL approximated as a whole-cell I_CaL current generated by cardiac-specific isoform Ca_v1.2 (i.e., percentage of Ca_v1.3 was set to 0). In this section, we report the results of our investigation of possible importance of Ca_v1.3 with respect to its interplay with CRU distribution. Because the contribution of Ca_v1.3 to the whole-cell I_CaL is presently unknown in rabbit SAN cells, we performed a full-scale parameter sensitivity analysis by varying contribution of Ca_v1.3 from 0 to 100% with an interval of 10%. Our simulations with Ca_v1.3 in I_CaL revealed the following.

(1) Substitution of Ca_v1.3 for Ca_v1.2 in the whole-cell I_CaL substantially shortened CL for both square lattice and uniformly random distributions of CRUs (Fig. 11 A). However, the total CL shortening with full substitution of Ca_v1.2 by Ca_v1.3 (i.e., 0 versus 100%) was smaller for the uniformly random CRUs, i.e., 25.9 versus 39.8% for square lattice CRUs (blue arrow versus orange arrow in Fig. 11 A).

(2) At any specific percentage of Ca_v1.3, the average CL was always notably shorter for the uniformly random CRU distribution. Thus, the CRU randomness caused an additional CL reduction on the top of CL reduction due to Ca_v1.3. This additional effect was substantially reduced as contribution of Ca_v1.3 increased. For example, with 0% of Ca_v1.3 in I_CaL, the CL shortening caused by CRU randomness was 31%, but only 15% with 100% of Ca_v1.3 in I_CaL (black arrow versus green arrow in Fig. 11 A).

(3) Mechanisms of CL reduction by Ca_v1.3 isoform and their interplay with effects of CRU randomness are shown in Fig. S4 and Fig. S5. The mechanisms are consistent with the concept of ignition (Lyashkov et al., 2018) that includes interactions of I_CaL, I_NCX, and CRUs (LCRs). For both square lattice and uniformly random distributions of CRUs, Ca_v1.3 substantially accelerated the timing and the rate of CRU recruitment that, in turn, accelerated I_NCX and diastolic depolarization. During diastolic depolarization, Ca_v1.3 (versus Ca_v1.2) not only causes stronger and earlier recruitment of CRUs to fire, but also generates a larger inward current due to its lower activation voltage of around −55 mV (versus −40 mV).

(4) Ca_v1.3 increases robustness of pacemaker function: substitution of only 20% of Ca_v1.2 by Ca_v1.3 in I_CaL was sufficient to revive spontaneous AP firing in the SAN cell model with decreased g_CaL and square lattice distribution of CRUs (Fig. S6) which was similar to the effect of CRU randomness described above (Fig. 7 C).

(5) We also tested whether the CRU randomness has any effect on βAR stimulation in the model with Ca_v1.3. This was tested in one model scenario with 50% of Ca_v1.3 in I_CaL. A similar substantial contribution (∼60%) has been demonstrated in mouse SAN cells (Christel et al., 2012) and we kept it substantial in this test in the absence of literature data for rabbit. Our simulations showed that while the average CL was substantially longer in the case of square lattice versus uniformly random distribution of CRUs (318.4 versus 260.8 ms; Fig. 11 B) in basal state firing, the CL converged to about the same level (247.8 versus 236.6 ms; Fig. 11 C) in the presence of βAR stimulation, clearly showing a stronger overall effect in case of the square lattice CRU distribution (22.2 versus 9.3% of CL reduction), i.e., qualitatively similar to our result with the I_CaL model with 100% Ca_v1.2 (Fig. 9).

Our previous results with local Ca control models (Maltsev et al., 2011; Maltsev et al., 2013; Stern et al., 2014; Maltsev et al., 2017b), the ignition theory of pacemaking (Lyashkov et al., 2018), and the results of the present study indicate the crucial importance of local, intimate interactions among CRUs, NCX, and LCCs, suggesting that clustering LCCs with RyRs can be important not only for normal excitation–contraction coupling in cardiac muscle cells, but also for cardiac pacemaking. Next, we performed a series of simulations to test this idea using a modification of our model, in which we artificially uncoupled Ca_v1.2 currents from CRUs. In these simulations, Ca_v1.2 currents were uniformly distributed over the cell membrane, while Ca_v1.3 currents remained coupled with CRUs. Our sensitivity analysis (similar to that shown in Fig. 11 A) revealed strong regulation of CL duration by both CRU randomness and Ca_v1.3 (Fig. S7). However, the robustness of AP firing in these cell models was substantially decreased. Indeed, while normal automaticity was present for the entire range of %Ca_v1.3 in the fully coupled channel models, the partially uncoupled models failed to generate spontaneous APs with %Ca_v1.3 being below 20% with random CRU distribution and below 30% with square lattice distribution (Fig. S7).

The focus of the present study is to investigate if disorder in CRU positions under cell membrane has any notable effect on pacemaker SAN function. We used an upgraded numerical SAN cell model featuring a Ca clock at the level of the local CRU network, coupled to the cell membrane electrophysiology equations. As the extent of disorder in the CRU positions increased, SAN function was assessed via its spontaneous AP firing rate or CL. Spatial disorder increased both the firing rate and robustness of pacemaker function in the basal state. The disorder decreased the CRU nearest neighbor distances and facilitated Ca release propagation via CICR. This leads to an earlier and stronger LCR signal that increased AP firing rate and also initiated AP firing in cells that could not fire APs with fully ordered square lattice positioning of CRU. The magnitude of the effect on the firing rate was substantial, quantitatively similar to that known for βAR stimulation. However, the boost of robustness bears a cost in the flexibility of pacemaker function. The range of AP firing rate modulation by βAR stimulation substantially narrowed as disorder and its attendant robustness increase. This happens because the disorder-facilitated release synchronization utilizes a majority of CRUs in the basal state and leaves a smaller fraction (i.e., a smaller reserve) of CRUs to be recruited in AP ignition during βAR stimulation.

An important result of the present study is a major upgrade of our previous CRU-based model of rabbit SAN cell (Maltsev et al., 2011; Maltsev et al., 2013; Fig.1; see Appendix for more details). Timing for CRU activation was previously defined phenomenologically via a fixed refractory period followed by a Poisson process, whereas termination was set simply to a mean value of spark duration. Both the refractory period and spark duration were taken directly from experimental measurements. The new model includes Ca-release activation and termination mechanisms linked SR Ca load as shown in experimental and theoretical studies (Zima et al., 2008; Maltsev and Lakatta, 2009; Imtiaz et al., 2010; Vinogradova et al., 2010; Zima et al., 2010; Laver et al., 2013; Stern et al., 2013; Maltsev et al., 2017a; Veron et al., 2021). These model enhancements are of crucial importance for this and future studies. Firstly, the new model reflects recent progress in our understanding of RyR function. Secondly, it describes mechanistically the coupling of the Ca clock to the membrane clock at the scale of the local CRU network. This is an important niche in the variety of numerical SAN cell models developed thus far. As described in details in the Materials and methods and Appendix, its submicron scale (e.g., Nivala et al., 2012 in ventricular myocytes) is positioned between simple common pool models (e.g., Kurata et al., 2002; Maltsev and Lakatta, 2009; Severi et al., 2012) and extremely complex RyR-based models (Stern et al., 2014; Maltsev et al., 2017b).

The models at each level of signal integration are important for understanding the big picture, i.e., how the phenomenon of heartbeat emerges, bridging the gap between scales (Clancy and Santana, 2020; Weiss and Qu, 2020) in the spirit of multiscale modeling (Qu et al., 2011). The CRU-based modeling of SAN cells has become especially helpful for future research with respect to recent discoveries of importance of local Ca signaling for SAN tissue function (Bychkov et al., 2020). It was shown that Ca signals are markedly heterogeneous in space, amplitude, frequency, and phase among cells comprising an HCN4⁺/CX43⁻ cell meshwork of SAN, and synchronized APs emerge from heterogeneous subcellular subthreshold Ca signals (modelled here; see Fig. 8). Cell heterogeneity and biological noise are key determinants of robust cardiac pacemaking (Guarina et al., 2022; Maltsev et al., 2022).

While a series of insightful multicellular models of SAN tissue have been developed (see for example Oren and Clancy, 2010; Inada et al., 2014; Gratz et al., 2018; Li et al., 2018), local Ca signaling has not been numerically studied at the tissue level yet. Our new CRU-based model has relatively low computational cost and it generates LCRs observed experimentally both in isolated cells and in intact SAN and, thus, it can be used as a functional unit in future multicellular modeling to investigate the role of local Ca signaling at the level of SAN tissue at the frontier of heart pacemaker research (Clancy and Santana, 2020; Weiss and Qu, 2020). 1 s of simulation of our model requires about 12 min of computation time via one thread of Intel Xeon W-2145 CPU @3.7 GHz. Modern video cards can run several thousands of processing units. For example, TITAN RTX features 4608 CUDA cores running at 1.77 GHz. If a single CUDA core can be programmed to run the present single cell model, then running thousands of such models in parallel would simulate a respective tissue model comprised of these many-cell models in reasonable time, like for the single-cell model presented here. New theoretical insights into SAN tissue function can be, in fact, achieved in tissue models comprised of as low as 49 cells (7 × 7 grid; Gratz et al., 2018). Our recent investigation of GPU-based model of SAN tissue included 625 cells (25 × 25 grid; Maltsev et al., 2022).

Empty spaces and clusters are an intrinsic feature of random spatial distributions known as Poisson clumping. Such emerging CRU clumps are clearly seen by eye in our examples of 2-D representations of CRU networks in Fig. 2 (left), and the clustering effect is quantitatively manifested by a broader distribution of nearest-neighbor distances with notably shorter averages (right). The disorder in CRU positions creates shortcuts, making it easier for released Ca to reach a neighboring CRU via CICR propagation and thus promoting CRU recruitment and synchronization. Once the Ca release becomes more synchronized, LCR sizes increase, the amplitude of LCR net diastolic signal also increases, and, very importantly, the timing of net LCR signal simultaneously shortens, as we previously demonstrated (Maltsev et al., 2011). Further effect of LCRs on AP firing rate is executed via NCX and I_CaL to accelerate diastolic depolarization (Figs. 4, 5, and 10) as postulated in the coupled-clock theory (Maltsev and Lakatta, 2009; Lakatta et al., 2010) and more recent ignition theory (Lyashkov et al., 2018), in line with previous numerical studies in a CRU-based SAN cell model (Maltsev et al., 2013) and RyR-based model (Stern et al., 2014).

The Ca release synchronization mechanism via local CRU recruitment is impacted by many factors including, for example, I_spark, i.e., Ca release flux of a single CRU which we studied previously (Maltsev et al., 2011). It is determined by SR Ca refilling kinetics, driven by PKA-dependent phosphorylation of Ca cycling proteins (Vinogradova et al., 2006; Lakatta et al., 2010). Mechanistically speaking, Ca clock ticks during diastole when SR refills to a certain threshold level so that Ca current via a single RyR channel becomes big enough to recruit its neighboring RyRs within CRU to generate a spark (Zima et al., 2010; Veron et al., 2021) and then I_spark becomes large enough to generate an LCR, i.e., a series of propagating sparks that depends on the local CRU distribution. Thus, the recruitment phase becomes delayed if the nearest-neighbor distances are larger (such as in square lattice arrangement) and require a larger SR Ca loading commensurate with larger I_spark to begin with. And vice versa, the recruitment starts earlier in the cycle if some of nearest-neighbor distances are shorter (such as in uniformly random arrangement).

With larger nearest-neighbor distances in the square lattice setting, many CRUs never fire during diastole in basal state beating, but become recruited during βAR stimulation (Fig. 10). For example, in basal state on average about 44% of CRUs fired at −35 mV (i.e., the end of diastolic depolarization) in the uniformly random model, but only 29% in the lattice model; however, during βAR stimulation about 53% of CRUs fired in either case, i.e., CRUs became equally well recruited independent of the model. Thus, the slower recruitment and the nonrecruited CRUs in the basal firing represent chronotropic reserve mechanisms that are utilized during βAR stimulation. The reserve is obviously limited by the number of CRUs, and if more CRUs are recruited to fire in each diastole in basal state, then the number of CRUs in reserve shrinks (and vice versa); that is why βAR stimulation effect is much smaller in the uniformly random model versus lattice model. This result is evidence of possible functional importance of the chronotropic reserve linked to CRU recruitment: we have previously demonstrated that CRU recruitment stabilizes diastolic I_NCX amplitude that explained, for example, a paradoxical effect of partial knockout of NCX in mice to reduce chronotropic reserve with no effect on the basal rate (Maltsev et al., 2013).

Thus far, CRU distribution in SAN cells has not been systematically studied. However, available data indicate that locations of CRUs in SAN cells do not form a perfect grid, but exhibit a notable degree of randomness (Lyashkov et al., 2007; Stern et al., 2014). Different SAN cells may have different CRU arrangement: central SAN cells have a substantial degree of disorder in their RyR cluster positions, whereas peripheral cells feature striations and more organized RyR clustering (Rigg et al., 2000; Musa et al., 2002). Our result may indicate that the central cells are capable of robustly generating high frequency signals, but unlikely to have a substantial CRU recruitment reserve to utilize during βAR stimulation (they may still have other mechanisms to increase their rate). On the contrary, the cells with more organized RyR clustering (wherever they are) may have a larger CRU recruitment reserve. Thus, the difference in spatial geometry of the CRU network in different parts of the SAN and, hence, the associated difference in βAR stimulation may contribute to the leading pacemaker site shifts during βAR stimulation (Brennan et al., 2020).

The real CRU distribution is neither uniformly random, nor perfectly spaced. By adjusting the CRU positioning, SAN cells can reach a balance of robustness and flexibility. This balance is required (and dictated) for each cell by its specific location and functional role within the cellular network of SAN tissue. This could be important for regulation of βAR stimulation response among individual cells within SAN tissue for its optimal integrated chronotropic response (Brennan et al., 2020; Yuan et al., 2020; Kim et al., 2021). Here, we show that the CL decrease during βAR response depends linearly on the CL in the basal state. This finding is line with recent experimental results that βAR stimulation synchronizes a broad spectrum of AP firing rates in SAN cells toward a higher population average (Kim et al., 2021).

Another possible importance of the disorder-facilitated AP firing could be regulation of cell dormancy, a recently discovered phenomenon manifested by the absence of automaticity of SAN cells that can be acquired via βAR stimulation (Kim et al., 2018; Tsutsui et al., 2018; Tsutsui et al., 2021). One important factor of cell dormancy is a lower density of I_CaL (Tsutsui et al., 2021). Here, we show that by rearranging its CRU positions, a SAN cell with a low I_CaL density can switch its functional state between being dormant or AP firing, as illustrated in Fig. 7. Nonfiring cells have also been recently found in intact, fully functional SAN (Bychkov et al., 2020; Fenske et al., 2020); their numbers substantially varied in different chronotropic states of the SAN (Fenske et al., 2020). A proposed new pacemaker mechanism involves synchronized cardiac impulses emerging from heterogeneous local Ca signals within and among cells of pacemaker tissue, including subthreshold signals in nonfiring cells (Bychkov et al., 2020). Here, we show that subthreshold V_m oscillations in case of uniformly random distribution are much more powerful and occurred at a faster frequency versus those in the square lattice model (Fig. 8). Thus, the spatial CRU distribution can regulate the capacity of cells to generate synchrony of impulses via subthreshold signaling.

Another interesting result of our sensitivity analysis is that the model with the perfect square lattice CRU locations (before it failed) was capable of generating AP firing with much longer CL (i.e., at a low rate of excitations), whereas in uniformly random scenario the model failed at a shorter CL (i.e., cannot reach lower rates; Fig. 7, B and C). Thus, pacemaker cells can harness this mechanism linked to CRU distribution to reach lower rates in addition to other known bradycardic mechanisms such as shift in voltage activation of I_f (DiFrancesco and Tromba, 1988), I_K,Ach, and protein dephosphorylation (decreasing clock coupling; Lyashkov et al., 2009). By having more regular spacing, SAN cells can generate APs at slower rates, without any additional specific Ca- or voltage-dependent mechanisms, which may be more threshold sensitive.

Pathological conditions and aging are usually associated with disorder in molecular positions, interactions, and functions. Here, we demonstrate that excessive disorder in CRU positions within SAN cells decreases fight-or-flight response while shrinking the range of lower AP rates (Figs. 7 B and 9), and hence can explain, in part, the limited range of heart rates associated with age and in disease. On the other hand, the mechanism of disorder-facilitated Ca release propagation can act together with increased sympathetic tone to compensate the age-intrinsic heart rate range decline with age (Tsutsui et al., 2016). With respect to atrial and ventricular cells, excessive disorder in CRU positions (in cardiac disease) is expected to facilitate Ca release propagation, i.e., waves formation, increasing the risk of life-threatening arrhythmia (Ter Keurs and Boyden, 2007).

While noise is a broad term that is usually associated with undesirable disturbances or fluctuations, many biological systems harness disorder to function. Randomness creates opportunities to exceed a threshold that is above the mean, analogous to the way a quantum particle can tunnel across a barrier while a classical deterministic particle cannot. It can improve signal transmission or detection, e.g., via stochastic resonance (McDonnell and Abbott, 2009). Random parameter heterogeneity among oscillators can consistently rescue the system from losing synchrony (Zhang et al., 2021). Randomness is critical for cardiac muscle cell function. The local control theory developed by Michael Stern in 1992 (Stern, 1992) predicted the Ca sparks (found later experimentally in Cheng et al., 1993) and explained smooth regulation of excitation–contraction coupling in cardiac muscle via statistics of success and failure of a CRU to generate a spark when LCCs open. We have recently shown that statistical physics approach is also helpful to understand spontaneous Ca spark activation and termination (via Ising formalism; Maltsev et al., 2017a; Maltsev et al., 2019; Veron et al., 2021). In the present study, we show that disorder could be also important for cardiac pacemaker function: disorder in CRU locations determines statistics of success and failures for a firing CRU to recruit to fire its neighboring CRU, observed as propagating LCRs that are critical for SAN cell function (Lakatta et al., 2010). Thus, disorder in CRU positions facilitates functional order in terms of LCR emergence via self-organization by means of positive feedback provided by CICR, culminating in higher AP firing rates, whereas order in CRU positions is associated with individual stochastic sparks, i.e., functional disorder for a major part of the diastolic depolarization duration, culminating in lower AP firing rates. A broader interpretation of our results is that disorder in a network featuring diffusion-reaction interactions can facilitate excitation propagation that may be applicable to RyR arrangement within a CRU to generate a spark (down scale) or cell-to-cell interactions in SAN tissue to generate cardiac impulse (up scale).

According to a modern ignition theory (Lyashkov et al., 2018), the diastolic depolarization is realized by positive feedback mechanisms among CRUs (i.e., LCRs), I_NCX, and I_CaL via their intertwined Ca and voltage dependencies. Ca currents play a key role in generating LCRs during diastolic depolarization via CICR (Huser et al., 2000; Bogdanov et al., 2001; Chen et al., 2009; Torrente et al., 2016). Therefore, in addition to I_CaL density, I_CaL voltage dependence (especially within the range of diastolic depolarization) must be important. Studies in mice demonstrated specific importance of the Ca_v1.3 LCC isoform with a lower voltage activation threshold (versus cardiac isoform Ca_v1.2) for pacemaker function (Mesirca et al., 2015). An indication of Ca_v1.3 general importance (rather than just specific to mouse) is that the loss of Ca_v1.3 function is associated with a human channelopathy linked to bradycardia (Baig et al., 2011). Because specific contribution of Ca_v1.3 in rabbit SAN cells is presently unknown and Ca_v1.3 expression shows a substantial cell-to-cell variability (recently shown in mouse SAN cells via immunolabeling technique; Louradour et al., 2022), we performed a full-scale sensitivity analysis varying Ca_v1.3 percentage from 0 to 100%.

Our simulations show that inclusion of Ca_v1.3 in I_CaL can provide a stronger and earlier recruitment of CRUs to fire during diastolic depolarization (Figs. S4 and S5) supporting and providing further insights to previous studies which proposed that Ca_v1.3 regulates RyR-dependent Ca release during SAN pacemaker activity (Christel et al., 2012; Torrente et al., 2016; Louradour et al., 2022). We also found interesting interplay between Ca_v1.3 and CRU randomness (Fig. 11 A; see details in Results). Their effects are complementary. CRU randomness generates additional CL shortening on top of the Ca_v1.3 effect (Fig. 11 A, green arrow), and Ca_v1.3 can also generate an additional CL shortening on the top CRU randomness (Fig. 11 A, blue arrow). On the other hand, both factors compete for CRU recruitment and accelerate AP firing via the same ignition process. Thus, each factor decreases the CL modulatory range of the other: the presence of randomness in CRU locations decreases the effect of Ca_v1.3, and vice versa, expression of Ca_v1.3 decreases the effect of the CRU randomness.

Increasing relative contribution of Ca_v1.3 boosts robustness of pacemaker function. For example, replacing only 20% of Ca_v1.2 by Ca_v1.3 in I_CaL can revive normal automaticity in a dormant cell with weak CICR among CRUs located in a square grid (Fig. S6). On the other hand, strong expression of Ca_v1.3 (and concomitant high basal rate) decreases the effect of βAR stimulation (compare Fig. 11, B and C, with Fig. 9 B). SAN cells with higher basal rates showed smaller responses to βAR stimulation also in experimental studies (Kim et al., 2021). The inhibitory effect of CRU randomness on βAR stimulation is preserved in cells with high expression of Ca_v1.3. For example, our model with 50% of Ca_v1.3 in I_CaL exhibited a much smaller effect of βAR stimulation in case of unfirmly random CRU distribution versus square lattice distribution (9.3 versus 22.2% of CL reduction).

Our simulations also showed that cell models with artificially uncoupled Ca_v1.2 currents from CRUs exhibited a larger modulatory range of Ca_v1.3 but less robust peacemaking (Fig. S7), pointing to functional importance of local crosstalk of CRUs with both Ca_v1.2 and Ca_v1.3 channels, despite Ca_v1.2 having a higher activation voltage threshold. Thus, clustering LCCs with RyRs can be of general importance, not only for normal excitation–contraction coupling in cardiac muscle, but also for cardiac pacemaking. Uncoupling LCC from RyR in pathological conditions and aging can deteriorate cardiac pacemaker function.

In our previous study, we showed that a CRU network lacking release propagation can acquire release propagation capability by introducing a subset of smaller “bridging” CRUs that create propagation shortcuts and allow sparks to jump from one firing CRU to its neighbor (Stern et al., 2014). In the present study, we show that bridging of CRU network is not required to achieve release propagation: the intrinsic disorder in CRU positions can naturally create the bridges and propagation shortcuts. On the other hand, our study is of a reductionist type focused on the disorder contribution, whereas the realistic CRU distribution in SAN cells has more complex, hierarchical structure that includes CRUs of various sizes (Stern et al., 2014). Thus, future studies will clarify the role of disorder in the more realistic settings with different CRU sizes and more precise CRU locations within the cell measured by super-resolution microscopy in 3-D (see pilot studies: Maltsev et al., 2016; Greiser et al., 2020) rather than in 2-D tangential sections measured by confocal microscopy (Stern et al., 2014). Furthermore, RyR distribution is dynamic (Asghari et al., 2020) and spacing between CRUs becomes shortened in failing hearts (Chen-Izu et al., 2007). Future studies on the cellular and molecular mechanisms regulating CRU distribution dynamics within cardiac cells will clarify how CRU order/disorder contributes to cell physiological and pathological function. While our numerical simulations clearly show a notable effect of disorder in CRU positions on Ca release synchronization and spontaneous AP firing, the theoretical mechanisms of this synchronization merit further studies. Based on our numerical study, one can envision that synchronization is happening as a critical phenomenon, with the criticality depending on model parameters, including spatial randomness of the CRUs. Possible approaches to study such systems with criticalities include a percolation phase transition or Ising formalism, similar to what has been recently suggested for RyR interactions via CICR within CRUs (Maltsev et al., 2017a; Maltsev et al., 2019). While we tested possible functional effects of Ca_v1.3, the exact contribution of Ca_v1.3 into total I_CaL and biophysical properties of Ca_v1.3 current in rabbit SAN cells are presently unknown and require further elucidation.

The rate of Ca release propagation is an important feature of both normal and abnormal Ca release signals in cardiac cells. Using numerical modeling, here we show that disorder in CRU locations increases the synchronization of Ca release in SAN pacemaker cells. This impacts on their pacemaker function via NCX current accelerating diastolic depolarization. While the disorder increases the rate and robustness of spontaneous AP firing, it simultaneously decreases βAR stimulation effect and the low range of lower rates. Our simulations also showed that Ca_v1.3 strongly recruits CRUs to fire during diastolic depolarization, thereby increasing AP firing rate and complementing effects of CRU distribution. Thus, order/disorder in CRU locations together with Ca_v1.3 expression regulates CRU recruitment and synchronization to fire during diastolic depolarization and could be harnessed by pacemaker cells to regulate their function. Excessive CRU disorder and/or overexpression of Ca_v1.3 boost pacemaker function but can limit heart rate range that may contribute to heart rate range decline with age and in disease.

In the present study, we performed a major update of our previous CRU-based numerical model of a central rabbit SAN cell (Maltsev et al., 2011; Maltsev et al., 2013). The model formulations for cell membrane currents are adopted from the Maltsev and Lakatta (2009) model that, in turn, stems from Kurata et al. (2002). The LCR in the model is approximated at the scale of an individual CRU that represents a cluster of Ca release channels (RyRs) embedded in the JSR, i.e., a Ca store located in close proximity to cell surface membrane, with only 20 nm of separation via a dyadic space. Each JSR is diffusely linked to the network of FSR that uptakes Ca from cytosol via SERCA pumping. Individual release channels are not modeled here, but we translate recent findings of RyR studies to the CRU level to introduce respective spark activation and termination mechanisms. All CRUs are identical and located in an equally spaced square grid under the cell membrane. The model allows each CRU position to vary around its original square lattice position to generate a variety of intermediate distributions with various degrees of disorder from perfect square lattice to uniformly random. The specific aims of the present study were to update the model and to investigate how the disorder in CRU positions influences SAN cell function.

We model a small SAN cell shaped as a cylinder of 53.28 μm in length and 6.876 μm in diameter. The cell membrane electrical capacitance of 19.8 pF is similar to that of 20 pF in the Zhang et al. (2000) model of a central SAN cell. Details of local Ca dynamics under the cell membrane are simulated on a submicron resolution (120 nm) square grid that divides the cell membrane and the submembrane space (dubbed subspace) into respective membrane patches and subspace voxels. Locations within the grid are defined by coordinates in the respective plane of the cell surface cylinder: along the cylinder (x axis) and around the cell cross-section (y axis). Our cell partitioning and respective voxel sizes to reproduce LCRs observed experimentally are schematically illustrated in the main text (Fig. 1). To avoid special considerations at the cell borders, the ends of the cylinder are connected (yielding a torus). Our cell compartments and voxel structure are essentially similar to that in the Stern et al. (2014) model that approximates Ca dynamics in three dimensions at the single RyR scale. However, we limited the cell partition to only three nested layers of voxels of substantially different scales reflecting respective essential Ca cycling components and processes happening at these scales (described below).

The first, most detailed level of approximation of local Ca dynamics is via 79,920 voxels (444 in x and 180 in y) of very small size 120 × 120 × nm under the entire cell membrane including the dyadic space (or cleft space) separating CRUs and the cell membrane. This thin layer of voxels describes local Ca release from individual CRUs, CRU-to-CRU interactions via Ca diffusion, and CRU interactions with cell membrane (including I_CaL, I_NCX, I_CaT, and I_bCa). Ca currents were computed for each membrane patch (120 × 120 nm) to generate local Ca fluxes contributing to local Ca dynamics.

The next approximation level of Ca dynamics is a deeper layer of voxels that have a larger size of 360 × 360 × 800 nm (ΔxΔyΔr) that includes the scale of JSR depth (60 nm). Each ring voxel has its cytoplasmic part and FSR part and some ring voxels are diffusively connected to JSRs (Fig. 1). While it appears that the geometric scale of ring voxel is of an order of magnitude larger than the JSR size, the actual FSR volume within each ring voxel is comparable with the JSR volume (described below in detail). Thus, this level of voxels describes the local dynamics of Ca transfer from FSR to JSR as well as enhanced local Ca pumping and diffusion fluxes due to close proximity to Ca release and Ca influx in its neighboring submembrane voxels.

In contrast to the Stern et al. (2014) model, the rest of the cell in our model does not have further geometric partitioning and it is lumped to one compartment, the “core,” where local Ca dynamics is less important. The core also has its cytosolic and FSR parts (as in ring voxels). Thus, it describes the bulk Ca uptake from cytosol to FSR and further accumulation, diffusion, and redistribution of the pumped and released Ca within the cell interior.

We place JSRs inside the respective ring voxels just below their outer side facing the cell membrane (Fig. 1). JSR volume (∼7.8 attoliter) is comparable with the volume of a ring voxel (∼91 attoliter) and the same volume cannot be occupied twice by different cell compartments. Therefore, volumes of the ring voxels are kept the same in the model by their extending into the core by the exact volume that the JSR occupies at their outer side. Thus, the actual core volume was calculated as the volume of the cylinder core decreased by the volume of all JSR volumes. We simulated different degrees of disorder of JSR positions by a random number generator within a Gaussian distribution along x and y (centered at the square lattice vertices) with a given SD that was the same for x and y directions. The model scenario with uniformly random positioning of CRU centers was also generated by a random number generator but with equal probability to occur in any submembrane voxel. JSR overlaps are excluded, i.e., any two JSRs cannot occupy the same cell volume.

As mentioned above, each ring voxel and the core are further partitioned into cytosol and FSR fractions. We modeled FSR as homogeneously distributed network within the cytosol. The cytosol fraction was set to 0.46 and FSR fraction to 0.035 (Stern et al., 2014). The remainder presumably contains myofilaments, mitochondria, nucleus, and other organelles. In turn, each FSR portion has capability to pump Ca locally from the respective cytosol portion of the same voxel, simulating local SERCA function. Because the submembrane voxels are extremely thin, only 20 nm depth, their contribution to Ca pumping and intra-SR diffusion are negligible and not modelled.

However, because FCC is calculated only once before the model run and does not carry any additional computational burden during actual simulations of Ca dynamics, we always used here the full approximation for the diffusion, i.e., more precise exponential decay, rather than a linear change over Δt. An advantage of this approach is that the model features more stable behavior in case we want to vary cell geometry, cell compartments, voxel sizes, or integration time (Δt).

Note 2: We have only a fraction of cell volume occupied with cytosol (or FSR). However, the same fraction will be for v and s in τ formulations and it cancels. The volume ratios v₂/(v₁ + v₂) and v₁/v₂ remain also unchanged because the fraction factor also cancels. Thus, all above formulations with formal geometric volumes are also valid for fractional volumes, assuming that the fraction of cytosol (or FSR) is evenly distributed within the volume.

The JSRs are located in close proximity to the cell membrane separated only by the layer of submembrane voxels of 20 nm depth, representing the dyadic space. Thus, each CRU releases Ca (described below) into its neighboring subspace voxels (occupying the dyadic space). Each CRU is refilled with Ca locally from FSR network via a fixed diffusional resistance. In the previous common pool models (Kurata et al., 2002; Maltsev and Lakatta, 2009), diffusion between JSR and FSR was described by a simple exponential transfer process with a fixed time constant (τ_tr = 40 ms). Here, we want to have a similar Ca transfer rate for each JSR, but locally. The τ_tr value in common pool models is for the whole-cell FSR volume that is substantially larger than JSR volume. Here, in the local model, each JSR is linked diffusively to FSR. Depending on its position, JSR can be connected to one ring voxel, two voxels, or four voxels. To get the respective share of the Ca flux, we split and distribute the Ca diffusion flow into nine elementary surface areas (120 × 120 nm) of the JSR (360 × 360 nm), each of which is connected to the respective ring voxel and transfer its Ca share with τ = 9 · 40 ms · v_FSR/(v_FSR+v_JSR) = 104.58 ms.

Note: Because the FSR fraction of 0.035 within cell volume is rather small, each relatively large ring voxel of 360 × 360 × 800 nm (∼91 attoliter) contains only a small FSR volume v_FSR = 3.186 attoliter. This is comparable (and even smaller) than the JSR volume in the model v_JSR = 7.776 attoliter.

It is important to note that our new model of spark initiation does not reflect an intrinsic time-dependent refractory process; rather, we implement here the idea that spark can occur only when SR gets refilled to a critical level that the current amplitude of individual RyR can initiate regenerative CICR among neighboring RyRs (Zima et al., 2010; Stern et al., 2013; Veron et al., 2021). Thus, the implementation of the new spark activation mechanism controlled by SR Ca refilling represents a major advance of our model, because the spark activation timing is now predicted by the model. Of note, our previous CRU-based SAN cell models (Maltsev et al., 2011; Maltsev et al., 2013) implemented SR Ca refiling contribution phenomenologically via a fixed parameter, the restitution time that was directly taken from experimental measurements.

We also introduced a new spark termination mechanism that is based on the current knowledge in this research area (Laver et al., 2013; Stern et al., 2013; Maltsev et al., 2017a), i.e., a Ca spark is generated via CICR among individual RyRs within a CRU and it sharply terminates due to induction decay (Laver et al., 2013) or a phase transition (similar to that known in Ising model; Maltsev et al., 2017a) when I_RyR(t), current via single RyR, becomes too small (due to JSR Ca depletion) to further support the CICR. The specific value of the critical current I_{spark_termination} is defined by the RyR interactions and beyond the capability of our CRU-based model. But the time when CICR wanes to the critical point is reflected by the amplitude of I_spark(t) being comparable with I_RyR(t) at a given JSR load, so that only one or a very few RyRs remain open at the termination time point when spark decays. Based on this logic, we tested a wide range of I_{spark_termination} to generate sparks of various durations and found that spark termination is well described with I_{spark_termination} set to a 0.175 pA. In our previous numerical model simulations and Ising theory of spark (Maltsev et al., 2017a), spark terminates when SR level depletes to a critical level of about 0.1 to 0.15 mM (for RyR clusters from 9 × 9 to 13 × 13). Unitary RyR current I_RyR becomes very small at these SR levels, namely within the range of 0.035–0.0525 pA, respectively, assuming I_RyR at 1 mM of SR Ca to be I_RyR,1mM = 0.35 pA (as in Stern et al., 2013). Thus, our chosen critical spark amplitude of 0.175 pA is reasonable for spark termination time point, reflecting only 2 or 3 remaining open RyR channels (of 144 total in our SAN cell model). As the SR continues to be further depleted of Ca, these few open channels cannot support any longer CICR within the CRU, and spark undergoes a sharp termination phase transition that is also in line with our numerical simulations of spark dynamics (Figs. 1 and 3 in Maltsev et al., 2017a). It is important to note that the spark termination mechanism does not include any intrinsic time-dependent inactivation, but simply reflects the drop of local Ca gradient over the JSR to the critical point that cannot sustain CICR among RyRs within the CRU (Laver et al., 2013; Stern et al., 2013; Maltsev et al., 2017a).

Cytosolic Ca is buffered by calmodulin (0.045 mM) throughout the cell: submembrane voxels, ring voxels, and the core. Each JSR features Ca buffering with calsequestrin (30 mM).

The SERCA uptake flux j_up is given by Eq. 1.

Electrophysiological formulations for cell membrane currents are adapted from Maltsev and Lakatta (2009). Major changes of the model include introduction of local Ca currents and modulation of local currents by local Ca. To introduce local currents and local modulation by Ca, the cell membrane is partitioned into small patches, with each patch facing its respective subspace voxel. We also omitted sustained inward current I_st and background Na current I_bNa because thus far the molecular identities for these currents have not been found and these currents are likely produced by NCX or other currents with known molecular identity (Lakatta et al., 2010). We adopted I_NCX density and Ca-dependent I_CaL inactivation for more realistic modulation by much higher local Ca concentrations under the cell membrane predicted by our local Ca control models in the present study (Fig. S1) and previous models (Maltsev et al., 2011; Maltsev et al., 2013; Stern et al., 2014; Maltsev et al., 2017b) reaching >10 μM versus common pool models (Kurata et al., 2002; Maltsev and Lakatta, 2009) predicting [Ca] in subspace only in sub-μM range in a bulk “subspace” compartment during diastole and 1–2 μM during Ca transient peak.

Ca_o = 2: Extracellular Ca concentration; K_o = 5.4: Extracellular K concentration; K_i= 140: Intracellular K concentration; Na_o = 140: Extracellular Na concentration; Na_i = 10: Intracellular Na concentration.

E_Na = E_T ∙ ln(Na_o/Na_i): Equilibrium potential for Na; E_K = E_T ∙ ln(K_o/K_i): Equilibrium potential for K; E_Ks= E_T∙ ln[(K_o+ 0.12∙ Na_o)/(K_i+ 0.12 ∙ Na_i)]: Reversal potential of I_Ks. Where E_T is “RT/F” factor = 26.72655 mV at 37°C, E_CaL= 45: Apparent reversal potential of I_CaL; E_CaT= 45: Apparent reversal potential of I_CaT.

L-type Ca current (I_CaL)

Ca-dependent I_CaL inactivation is described by parameters K_mfCa and α_fCa. With K_mfCa of 0.35 μM in previous common pool models, our local Ca model of SAN cell did not work (not enough I_CaL was available to generate normal APs), because I_CaL inactivation was too sensitive to Ca and a major part of I_CaL was inactivated in diastole by diastolic LCRs whose amplitudes reach tens and hundreds of micromolars. Thus, we adopt Ca-dependent inactivation of I_CaL for more realistic modulation by much higher local Ca concentrations by setting K_mfCa to 30 μM.

We set the midpoint of I_CaL activation V_dL to −6.6 mV as in SAN cell models of Wilders et al. (1991) and Dokos et al. (1996). This value of V_dL is relatively high with respect to other SAN cell models and therefore was considered in our study to simulate I_CaL generated by Ca_v1.2, the cardiac isoform of LCC. It was used in all our simulations, except in the last two Results sections where we studied effects of inclusion of Ca_v1.3 into I_CaL. Both isoforms Ca_v1.2 and Ca_v1.3 are colocalized with RyRs (Christel et al., 2012), but Ca_v1.3 has a lower (more hyperpolarized) voltage activation threshold. Therefore, we simulated I_CaL generated by Ca_v1.3 using the same formulations above, but V_dL was set to −13.5 mV. The resultant V_dL shift increased d_L,∞ of Ca_v1.3 within the range of diastolic depolarization by a factor of approximately three, generating a much larger current and stronger recruitment of CRUs to fire during diastolic depolarization. We performed a full-scale sensitivity analysis varying percentage of Ca_v1.3 in g_CaL from 0 to 100%, with total conductance g_CaL of Ca_v1.2 and Ca_v1.3 remaining constant (gCa_v1.2 + gCa_v1.3 = g_CaL = 0.464 nS/pF = const).

T-type Ca current (I_CaT)

The whole-cell I_CaT was evenly distributed over cell membrane patches to generate respective homogeneous Ca influx. In each submembrane i-th voxel, I_CaT,i= I_CaT/N_voxels

Rapidly activating delayed rectifier K⁺current (I_Kr)

Slowly activating delayed rectifier K⁺current (I_Ks)

4-aminopyridine-sensitive currents (I_4AP= I_to+ I_sus)

Hyperpolarization-activated, “funny” current (I_f)

Na⁺-K⁺pump current (I_NaK)

The whole-cell I_bCa was evenly distributed over cell membrane to generate respective homogeneous Ca influx. In each submembrane i-th voxel, I_bCa,i= I_bCa/N_voxels.

I_NCX

Initial values for electrophysiology are in Table A1.

Initial values for Ca dynamics are in Table A2.

Ca_o = 2 mM: Extracellular [Ca]; K_o = 5.4 mM: Extracellular [K]; Na_o = 140 mM: Extracellular [Na]; K_i = 140 mM: Intracellular [K]; Na_i = 10 mM: Intracellular [Na].

E_CaL = 45 mV: Apparent reversal potential of I_CaL; g_CaL = 0.464 nS/pF: Conductance of I_CaL; K_mfCa = 0.03 mM: Dissociation constant of Ca-dependent I_CaL inactivation; beta_fCa = 60 mM⁻¹ · ms⁻¹: Ca association rate constant for I_CaL; alfa_fCa = 0.021 ms⁻¹: Ca dissociation rate constant for I_CaL; k_tau_fL = 0.5: Scaling factor for tau_fL used to tune the model; E_CaT = 45 mV: Apparent reversal potential of I_CaT, mV; g_CaT = 0.1832 nS/pF: Conductance of I_CaT; g_If = 0.105 nS/pF: Conductance of If; V_If,1/2 = −64 mV: half activation voltage of I_f; g_Kr = 0.05679781 nS/pF: Conductance of delayed rectifier K current rapid component; k_tau_I_Kr = 0.3: scaling factor for tau_paF and tau_paS used to tune the model; g_Ks = 0.0259 nS/pF: Conductance of delayed rectifier K current slow component; g_to = 0.252 nS/pF: Conductance of 4-aminopyridine sensitive transient K current; g_sus = 0.02 nS/pF: Conductance of 4-aminopyridine sensitive sustained K current; I_NaKmax = 1.44 pA/pF: Maximum Na⁺/K⁺ pump current; K_mKp = 1.4 mM: Half-maximal Ko for I_NaK; K_mNap = 14 mM: Half-maximal Na_i for I_NaK; g_bCa = 0.003 nS/pF: Conductance of background Ca current; k_NCX = 48.75 pA/pF: Maximum amplitude of I_NCX.

K_1ni = 395.3 mM: intracellular Na binding to first site on NCX; K_2ni = 2.289 mM: intracellular Na binding to second site on NCX; K_3ni = 26.44 mM: intracellular Na binding to third site on NCX; K_1no = 1628 mM: extracellular Na binding to first site on NCX; K_2no = 561.4 mM: extracellular Na binding to second site on NCX; K_3no = 4.663 mM: extracellular Na binding to third site on NCX; K_ci = 0.0207 mM: intracellular Ca binding to NCX transporter; K_co = 3.663 mM: extracellular Ca binding to NCX transporter; K_cni = 26.44 mM: intracellular Na and Ca simultaneous binding to NCX.

Q_ci = 0.1369: intracellular Ca occlusion reaction of NCX; Q_co = 0: extracellular Ca occlusion reaction of NCX; Q_n = 0.4315: Na occlusion reactions of NCX.

k_bCM = 0.542 ms⁻¹: Ca dissociation constant for calmodulin; k_fCM = 227.7 mM⁻¹· ms⁻¹: Ca association constant for calmodulin; k_bCQ = 0.445 ms⁻¹: Ca dissociation constant for calsequestrin; k_fCQ = 0.534 mM⁻¹· ms⁻¹: Ca association constant for calsequestrin; CQ_tot = 30 mM: Total calsequestrin concentration; CM_tot = 0.045 mM: Total calmodulin concentration.

K_mf = 0.000246 mM: the cytosolic side K_d of SR Ca pump; K_mr = 1.7 mM: the lumenal side K_d of SR Ca pump; H = 1.787: cooperativity of SR Ca pump; P_up = 0.014 mM/ms: Maximal rate of Ca uptake by SR Ca pump.

CRU_Casens = 0.00015 mM: sensitivity of Ca release to Ca_sub; CRU_ProbConst = 0.00027 ms⁻¹: CRU open probability rate at Ca_sub = CRU_Casens; CRU_ProbPower = 3: Cooperativity of CRU activation by Ca_sub; CaJSR_spark_activation = 0.3 mM: critical JSR Ca loading to generate a spark (CRU can open); Ispark_Termination = 0.175 pA: critical minimum Ispark triggering spark termination (CRU closes); Ispark_activation_tau_ms = 80 ms: Time constant of spark activation (a); Iryr_at_1 mM_CaJSR = 0.35 pA: unitary RyR current at 1 mM delta Ca; RyR_to_RyR_distance_um = 0.03 μm: RyR crystal grid size in JSR; JSR_depth_um = 0.06 μm: JSR depth; JSR_Xsize_um = 0.36 μm: JSR size in x; JSR_Ysize_um = 0.36 μm: JSR size in y; JSR_to_JSR_X_um = 1.44 μm: JSR crystal grid size in x; JSR_to_JSR_Y_um = 1.44 μm: JSR crystal grid size in y.

L_cell = 53.28 μm: Cell length; r_cell = 3.437747 μm: Cell radius; C_m = 19.80142 pF: membrane electrical capacitance of our cell model with 0.0172059397937 pF/μm2 specific membrane capacitance calculated from Kurata et al. (2002) model for its 32 pF cylinder cell of 70 μm length and 4 μm radius; 0.12 μm: the grid size; 0.12 μm: Submembrane voxel size in x and y; 0.02 μm: Submembrane voxel size depth; 0.36 μm: Ring voxel size in x; 0.36 μm: Ring voxel size in y; 0.8 μm: Ring voxel depth; 0.46: Fractional volume of cytosol; 0.035: Fractional volume of FSR.

D_cyt = 0.35 μm²/ms: Diffusion coefficient of free Ca in cytosol; D_FSR = 0.06 μm²/ms: Diffusion coefficient of free Ca in FSR.

For square lattice and uniform random distributions of CRUs, we performed sensitivity analysis for g_CaL from its basal value of 0.464 nS/pF (100%) down to 0.2552 nS/pF (55%) shown by magenta band in Fig. 7 A with a step of 0.0232 pA/pF (5%). The original data of this analysis in the form of V_m time series are given Fig. S2 for each CRU distribution.

Effect of βAR stimulation was modeled essentially as we previously reported (Maltsev and Lakatta, 2010) by increasing I_CaL, I_Kr, I_f, and Ca uptake rate by FSR via SERCA pumping. Specifically, whole-cell maximum I_CaL conductance g_CaL was increased by a factor of 1.75 from 0.464 to 0.812 nS/pF; whole-cell maximum I_Kr conductance g_Kr was increased by a factor of 1.5 from 0.05679781 to 0.085196715 nS/pF; the midpoint of I_f activation curve was shifted to more depolarized potential by 7.8 mV from −64 to −56.2 mV; and the maximum Ca uptake rate P_up was increased by a factor of 2 from 14 to 28 mM/s. The original data in the form of intervalograms are given (Fig. S3) for each CRU distribution in basal state and in during βAR stimulation.

The model code was written in Delphi Language (Delphi 10.4) and was computed with a fixed time tick of 0.0075 ms on a workstation running Windows 10 with Intel Xeon W-2145 CPU @3.7 GHz processor. The basic model code of the model is provided as a single Delphi file (SANC.dpr; Data S1) that can be freely used as a new mainframe to further investigate SAN function with respect to local Ca changes and Ca channel isoforms Ca_v1.2 and Ca_v1.3 interacting locally with RyRs.

David A. Eisner served as editor.

This research was supported in part by the Intramural Research Program of the National Institutes of Health, National Institute on Aging. A.V. Maltsev acknowledges the support of the Royal Society University Research Fellowship UF160569.

The authors declare no competing financial interests.

Author contributions: A.V. Maltsev, V.A. Maltsev, M.D. Stern: Conceptualization, Methodology, Validation.  A.V. Maltsev, V.A. Maltsev:  Software, Formal analysis, Investigation, Writing - Original Draft. A.V. Maltsev, M.D. Stern: Funding acquisition, Resources, Project administration.  V.A. Maltsev: Data Curation, Visualization.