Excitation–contraction coupling kinetics is dictated by the action potential rate of sinoatrial-nodal cells. These cells generate local Ca releases (LCRs) that activate Na/Ca exchanger current, which accelerates diastolic depolarization and determines the pace. LCRs are generated by clusters of ryanodine receptors, Ca release units (CRUs), residing in the sarcoplasmic reticulum. While CRU distribution exhibits substantial heterogeneity, its functional importance remains unknown. Using numerical modeling, here we show that with a square lattice distribution of CRUs, Ca-induced-Ca-release propagation during diastolic depolarization is insufficient for pacemaking within a broad range of realistic ICaL densities. Allowing each CRU to deviate randomly from its lattice position allows sparks to propagate, as observed experimentally. As disorder increases, the CRU distribution exhibits larger empty spaces and simultaneously CRU clusters, as in Poisson clumping. Propagating within the clusters, Ca release becomes synchronized, increasing action potential rate and reviving pacemaker function of dormant/nonfiring cells. However, cells with fully disordered CRU positions could not reach low firing rates and their β-adrenergic–receptor stimulation effect was substantially decreased. Inclusion of Cav1.3, a low-voltage activation L-type Ca channel isoform into ICaL, strongly increases recruitment of CRUs to fire during diastolic depolarization, increasing robustness of pacemaking and complementing effects of CRU distribution. Thus, order/disorder in CRU locations along with Cav1.3 expression regulates pacemaker function via synchronization of CRU firing. Excessive CRU disorder and/or overexpression of Cav1.3 boosts pacemaker function in the basal state, but limits the rate range, which may contribute to heart rate range decline with age and disease.

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 (INCX) 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 Ispark; 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 Cav1.3 LCC isoform for cardiac pacemaker function (Mesirca et al., 2015). While both Cav1.2 and Cav1.3 isoforms are colocalized with RyR clusters (i.e., CRUs; Christel et al., 2012), Cav1.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 Cav1.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 Ispark 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 Ispark 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) ICaL formulation was improved by including a contribution of Cav1.3 isoform to ICaL in addition to the cardiac-specific isoform Cav1.2. Both Cav1.2 and Cav1.3 isoforms are localized in CRUs (Christel et al., 2012), but have different voltage activation thresholds (see Appendix). Our simulations were initially performed with ICaL comprised of cardiac-specific isoform Cav1.2. Then, we investigated effects of inclusion of Cav1.3 isoform into ICaL (see the last two sections of the Results).

Online supplemental material

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 Cav1.3 accelerates AP firing rate by interacting with CRUs and INCX in the model with square lattice distribution of CRUs. Fig. S5 shows that Cav1.3 accelerates AP firing rate by interacting with CRUs and INCX in the model with uniformly random distribution of CRUs. Fig. S6 shows that Cav1.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.

The basal AP firing rate increases as CRU distribution becomes disordered

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% Cav1.2 in ICaL) 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 (R2 = 0.96; Fig. 3 A, inset), indicating importance of CRU interactions with their nearest neighbors (i.e., CICR) in this phenomenon.

Disorder in CRU positions accelerates CRU recruitment and Ca release synchronization in diastole

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 NCRU, 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.

Electrophysiological consequences of accelerated CRU recruitment: Role of INCX and ICaL

As a result of faster CRU recruitment, the INCX and ICaL were activated earlier in the cycle, accelerating diastolic depolarization and thereby shortening the CL (Fig. 4, C and D). In all scenarios, ICaL was activated concurrently with INCX 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 INCX at low voltages, below ICaL 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 (Vm) of −45 mV that is close to the AP ignition onset and ICaL 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).

Video 1.

Square lattice distribution of CRUs. Simulation of Ca dynamics in submembrane space ([Ca]sub) during 1 s. Note mainly individual sparks in diastole. [Ca]sub is coded from 0.15 μM (black) to 10 μM (saturation) by a color scheme shown at the bottom of Fig. 6; open CRUs are shown by white dots. Closed CRUs in refractory period are shown by blue dots; closed reactivated CRUs (available to fire) are shown by green dots. JSR Ca level is coded by respective shade (white, blue, or green) with a saturation level set at 0.3 mM. Simulation time and membrane potential (Vm) are shown in the top left corner.

Video 1.

Square lattice distribution of CRUs. Simulation of Ca dynamics in submembrane space ([Ca]sub) during 1 s. Note mainly individual sparks in diastole. [Ca]sub is coded from 0.15 μM (black) to 10 μM (saturation) by a color scheme shown at the bottom of Fig. 6; open CRUs are shown by white dots. Closed CRUs in refractory period are shown by blue dots; closed reactivated CRUs (available to fire) are shown by green dots. JSR Ca level is coded by respective shade (white, blue, or green) with a saturation level set at 0.3 mM. Simulation time and membrane potential (Vm) are shown in the top left corner.

Close modal
Video 2.

Perturbed square lattice distribution of CRUs with SD = 0.25 μm. Sizes of LCRs increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Video 2.

Perturbed square lattice distribution of CRUs with SD = 0.25 μm. Sizes of LCRs increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Close modal
Video 3.

Perturbed square lattice distribution of CRUs with SD = 0.5 μm. Sizes of LCRs further increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Video 3.

Perturbed square lattice distribution of CRUs with SD = 0.5 μm. Sizes of LCRs further increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Close modal
Video 4.

Perturbed square lattice distribution of CRUs with SD = 0.75 μm. Sizes of LCRs further increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Video 4.

Perturbed square lattice distribution of CRUs with SD = 0.75 μm. Sizes of LCRs further increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Close modal
Video 5.

Uniformly random distribution of CRUs. Sizes of LCRs further increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Video 5.

Uniformly random distribution of CRUs. Sizes of LCRs further increase during diastolic depolarization via propagating CICR. See Video 1 legend for color description and other details.

Close modal

Disorder in CRU positions increases robustness of basal automaticity

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 ICaL (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 ICaL densities at which spontaneous AP firing is impossible for the square lattice CRU distribution, but possible for uniformly random model. To match realistic ICaL density values (in pA/pF) and our model parameter gCaL (in nS/pF) determining maximum ICaL 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 ICaL 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 ICaL densities approximately in the middle. Please note that the graph in Fig. 7 A has a dual y-axis scale for ICaL and gCaL, respectively, bridging experimental and simulated data.

Then, we performed a wide-range sensitivity analysis for gCaL 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 ICaL values (magenta band in Fig. 7 A). The analysis revealed a wide margin of gCaL 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 gCaL 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.

Subthreshold signaling

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 ICaL (Tsutsui et al., 2021). Using the results of our sensitivity analysis, we chose a low gCaL of 0.2552 nS/pF, at which AP firing is impossible in both uniformly random and square lattice models and compared subthreshold Vm signals in the nonfiring cells in each case (Fig. 8 A). We found that subthreshold Vm 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 Vm signals computed for the time period after AP failure (Fig. 8 B).

Video 6.

Stochastic individual sparks generated by a dormant cell model with square lattice distribution of CRUs. The dormancy was achieved by decreasing ICaL conductance gCaL to 0.2552 nS/pF. Simulation time and membrane potential (Vm) are shown in the top left corner. See Video 1 legend for color description and other details.

Video 6.

Stochastic individual sparks generated by a dormant cell model with square lattice distribution of CRUs. The dormancy was achieved by decreasing ICaL conductance gCaL to 0.2552 nS/pF. Simulation time and membrane potential (Vm) are shown in the top left corner. See Video 1 legend for color description and other details.

Close modal
Video 7.

Subthreshold oscillatory LCR signals generated by a dormant cell model with uniformly random distribution of CRUs. The dormancy was achieved by decreasing ICaL conductance gCaL to 0.2552 nS/pF. Simulation time and membrane potential (Vm) are shown in the top left corner. See Video 1 legend for color description and other details.

Video 7.

Subthreshold oscillatory LCR signals generated by a dormant cell model with uniformly random distribution of CRUs. The dormancy was achieved by decreasing ICaL conductance gCaL to 0.2552 nS/pF. Simulation time and membrane potential (Vm) are shown in the top left corner. See Video 1 legend for color description and other details.

Close modal

Order in CRU position increases chronotropic reserve in the fight-or-flight reflex

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 (R2 = 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 NCRU and INCX 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 INCX. 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.

Interplay of Cav1.3 and CRU distribution

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

(1) Substitution of Cav1.3 for Cav1.2 in the whole-cell ICaL 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 Cav1.2 by Cav1.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 Cav1.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 Cav1.3. This additional effect was substantially reduced as contribution of Cav1.3 increased. For example, with 0% of Cav1.3 in ICaL, the CL shortening caused by CRU randomness was 31%, but only 15% with 100% of Cav1.3 in ICaL (black arrow versus green arrow in Fig. 11 A).

(3) Mechanisms of CL reduction by Cav1.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 ICaL, INCX, and CRUs (LCRs). For both square lattice and uniformly random distributions of CRUs, Cav1.3 substantially accelerated the timing and the rate of CRU recruitment that, in turn, accelerated INCX and diastolic depolarization. During diastolic depolarization, Cav1.3 (versus Cav1.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) Cav1.3 increases robustness of pacemaker function: substitution of only 20% of Cav1.2 by Cav1.3 in ICaL was sufficient to revive spontaneous AP firing in the SAN cell model with decreased gCaL 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 Cav1.3. This was tested in one model scenario with 50% of Cav1.3 in ICaL. 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 ICaL model with 100% Cav1.2 (Fig. 9).

Importance of L-type channel coupling to CRUs

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 Cav1.2 currents from CRUs. In these simulations, Cav1.2 currents were uniformly distributed over the cell membrane, while Cav1.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 Cav1.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 %Cav1.3 in the fully coupled channel models, the partially uncoupled models failed to generate spontaneous APs with %Cav1.3 being below 20% with random CRU distribution and below 30% with square lattice distribution (Fig. S7).

Results summary

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.

Upgraded CRU-based model of SAN cell

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).

Mechanisms of the CRU spatial disorder effect

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 ICaL 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, Ispark, 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 Ispark 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 Ispark 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 INCX 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).

Possible importance for normal function, pathological conditions, and aging

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 ICaL (Tsutsui et al., 2021). Here, we show that by rearranging its CRU positions, a SAN cell with a low ICaL 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 Vm 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 If (DiFrancesco and Tromba, 1988), IK,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).

A broader interpretation: Living systems harness disorder to function

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).

Importance of Cav1.2 and Cav1.3 isoforms

According to a modern ignition theory (Lyashkov et al., 2018), the diastolic depolarization is realized by positive feedback mechanisms among CRUs (i.e., LCRs), INCX, and ICaL 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 ICaL density, ICaL voltage dependence (especially within the range of diastolic depolarization) must be important. Studies in mice demonstrated specific importance of the Cav1.3 LCC isoform with a lower voltage activation threshold (versus cardiac isoform Cav1.2) for pacemaker function (Mesirca et al., 2015). An indication of Cav1.3 general importance (rather than just specific to mouse) is that the loss of Cav1.3 function is associated with a human channelopathy linked to bradycardia (Baig et al., 2011). Because specific contribution of Cav1.3 in rabbit SAN cells is presently unknown and Cav1.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 Cav1.3 percentage from 0 to 100%.

Our simulations show that inclusion of Cav1.3 in ICaL 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 Cav1.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 Cav1.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 Cav1.3 effect (Fig. 11 A, green arrow), and Cav1.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 Cav1.3, and vice versa, expression of Cav1.3 decreases the effect of the CRU randomness.

Increasing relative contribution of Cav1.3 boosts robustness of pacemaker function. For example, replacing only 20% of Cav1.2 by Cav1.3 in ICaL 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 Cav1.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 Cav1.3. For example, our model with 50% of Cav1.3 in ICaL 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 Cav1.2 currents from CRUs exhibited a larger modulatory range of Cav1.3 but less robust peacemaking (Fig. S7), pointing to functional importance of local crosstalk of CRUs with both Cav1.2 and Cav1.3 channels, despite Cav1.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.

Limitations and future studies

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 Cav1.3, the exact contribution of Cav1.3 into total ICaL and biophysical properties of Cav1.3 current in rabbit SAN cells are presently unknown and require further elucidation.

Conclusion

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 Cav1.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 Cav1.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 Cav1.3 boost pacemaker function but can limit heart rate range that may contribute to heart rate range decline with age and in disease.

1. General description of the model

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.

2. Cell geometry, compartments, voxels, and membrane patches

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).

2.1 Submembrane voxels

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 ICaL, INCX, ICaT, and IbCa). Ca currents were computed for each membrane patch (120 × 120 nm) to generate local Ca fluxes contributing to local Ca dynamics.

2.2 JSR level voxels (dubbed ring voxels)

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.

2.3 The core

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.

2.4 JSR

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.

3. Ca cycling

3.1 FSR

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.

3.2 SR Ca pump

The SERCA pump is present uniformly throughout the cell, transferring Ca from the cytosolic to the FSR compartment of each voxel (of ring and core) with Ca uptake rate given by the reversible Ca pump formulation adopted from Shannon et al. (2004).
$jup=PupVmax⋅(CacytKmf)H−Vmax⋅(CaFSRKmr)H1+(CacytKmf)H+(CaFSRKmr)H.$
(1)
where Pup = 0.014 mM/ms, Kmf = 0.000246 mM, Kmr = 1.7 mM, and H = 1.787.

3.3 Ca diffusion within and among cell compartments

Ca diffusion fluxes between voxels within and among cell compartments are approximated by the first Fick’s law:
$J=−D·Δ[Ca]/Δx.$
where D is a diffusion coefficient, and Δ[Ca]/Δx is Ca concentration gradients, i.e., Δ[Ca] is the concentration difference and Δx is the distance between the voxel centers.
The respective rate of change of [Ca] is defined as
$dCa/dt=J·s/v.$
where s is the diffusion area sharing by voxels and v is the receiving volume. Thus,
$dCa/dt=−D·(Δ[Ca]/Δx)·s/v.$
For any two diffusively interacting voxels with volumes v1 and v2, Ca dynamics is described by a set of differential equations:
$dCa1/dt=(Ca2−Ca1)/τ1,$
$dCa2/dt=(Ca1−Ca2)/τ2.$
where Ca1 and Ca2 are Ca concentrations in the respective voxels and
$τ1=v1·Δx/(D·s),$
$τ2=v2·Δx/(D·s).$
τ1 (or τ2, symmetrically) is the respective time constant of Ca1 change in time in a special case if the other compartment with volume v2 is substantially larger than v1 (i.e., v2v1) and therefore Ca2 remains approximately constant. In general case, the set is analytically solved to the respective exponential decays:
$Ca1(t)=(Ca1(0)−C∞)·exp(−tτ)+C∞,$
$Ca2(t)=(Ca2(0)−C∞)·exp(−tτ)+C∞.$
where C = (Ca1(0)·v1 + Ca2(0)·v2)/(v1 + v2) is equilibrium concentration (t = ∞) in both voxels defined by the matter conservation principle and τ = τ1 · τ2/(τ1 + τ2) is the common time constant of the exponential decay of the system to reach the equilibrium. The respective Ca change from its initial value in voxel v1 over time is given as follows:
$ΔCa1t=C∞−Ca10·1−exp−tτ.$
Then, by substituting C we get:
$ΔCa1(Δt)=[v2(v1+v2)]·(Ca2(0)−Ca1(0))·[1−exp(−Δtτ)],$
$ΔCa2(Δt)=−[v1(v1+v2)]·(Ca2(0)−Ca1(0))·[1−exp(−Δtτ)]=−ΔCa1(t)·v1v2.$
These formulations were used in all our computations of [Ca] changes for all neighboring voxels within and among cell compartments for the model integration for each time update Δt (during time tick or several time ticks for slower processes). In our computer algorithm, we calculated the fractional Ca change (FCC) before the model run and used it simply as a scaling factor to determine actual [Ca] change from the difference in [Ca] between any two diffusely interacting voxels at the beginning of each integration step (from t = 0 to t = Δt). Thus,
$FCC=[v2/(v1+v2)]·[1−exp(−Δtτ)]$
(calculated before model run),
$ΔCa1(Δt)=FCC·(Ca2(0)−Ca1(0))$
(calculated during model run), and
$ΔCa2(Δt)=−ΔCa1(t)·v1v2$
(calculated during model run).
In the case of identical voxels, i.e., within subspace and within ring (i.e., when v1= v2 and τ1= τ2) the formulations are simplified to
$FCC=0.5·[1−exp(−2·Δtτ1)].$
Note 1: If Δt/τ1 ≪ 1, then FCC can be approximated (e.g., via respective Taylor series) as
$FCC=Δtτ1=Δt·D·sv1·Δx.$
Further, if v1 can be described as v1= s · Δx, e.g., for diffusion along the cell length (axis x in our model), FCC can be further simplified to
$FCC=Δt·D/Δx2.$

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 v2/(v1 + v2) and v1/v2 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.

3.4 JSR and Ca diffusion between JSR and FSR

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 · vFSR/(vFSR+vJSR) = 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 vFSR = 3.186 attoliter. This is comparable (and even smaller) than the JSR volume in the model vJSR = 7.776 attoliter.

3.5 Spark activation mechanism

Each CRU can be either in open or closed state. The capability of a given CRU to open, i.e., to generate a Ca spark, is controlled by its JSR Ca loading. Experimental and theoretical studies showed that sparks cannot be generated with SR Ca loading less than a certain critical level of about 300 μm (Zima et al., 2010; Veron et al., 2021). This critical level of CajSR is implemented in our mechanism of spark activation by prohibiting CRU firing while SR Ca loading remains below 300 μM (CaJSR_spark_activation). When JSR is refilled with Ca above this level, it can open. The switch from close state to open state is probabilistic. The probability density for a given closed CRU to open is described by a power function of Ca concentration (Ca) in the dyadic space. The probability for the CRU to open during a short time interval, TimeTick, is given by
$P=ProbConst·(Ca/Casens)ProbPower·TimeTick.$
(2)
where Casens = 0.15 μM sensitivity of CRU to Ca, ProbConst = 0.00027 ms−1 is open probability rate at Ca = Casens, and ProbPower = 3 defines the cooperativity of CRU activation by cytosolic Ca. Each time tick, our computer algorithm tries to activate a closed CRU by generating a random number within (0, 1). If this number is less than P, then the CRU opens. When CRU opens its Ca current amplitude, Ispark, is defined by spark activation kinetics a(t), RyR unitary current (IRyR,1mM, the current via a single RyR at 1 mM of Ca gradient), the number of RyRs residing in the JSR (NRyR), and concentration difference between inside and outside JSR:
$Ispark=a(t)·IRyR,1mM·NRyR·(Ca−CaJSR).$
(3)
NRyR is defined based on the surface area of JSR assuming a crystal-like structure for RyR positions separated by 30 nm. Thus, for our JSR xy area of 360 × 360 nm, we obtain 12 × 12 RyRs, i.e., NRyR = 144. IRyR,1mM is to set 0.35 pA (Stern et al., 2014). Spark activation is described as a single exponential time-dependent process to tune spark rise time to about 5 ms (Fig. S1) close to that reported in the literature (from 4 to 8 ms).
$a(t)=1−exp(−t/τCRUactivation).$
(4)

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.

3.6 Spark termination mechanism

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 IRyR(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 Ispark_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 Ispark(t) being comparable with IRyR(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 Ispark_termination to generate sparks of various durations and found that spark termination is well described with Ispark_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 IRyR becomes very small at these SR levels, namely within the range of 0.035–0.0525 pA, respectively, assuming IRyR at 1 mM of SR Ca to be IRyR,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).

3.7 Ca buffering

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).

3.8 Summary of equations of local Ca dynamics

In a submembrane voxel:
$∂Cacyt,sub∂t=Dcyt∇2Cacyt,sub−jcyt_sub_to_ring+Ispark/Nvoxels_in_dyad−iCa/2·F·vcyt,sub−JCM.$
(5)
Ispark in voxels outside dyadic space is absent. Each CRU releases Ca (given by Ispark) into its dyadic space. Ispark is evenly distributed among submembrane voxels of the dyadic space (Nvoxels_in_dyad= 9). Ispark is given in Eq. 3 and it is positive, i.e., increasing [Ca] in the submembrane voxel. iCa is the sum of local Ca transmembrane currents (described in detail below in Electrophysiology) via the membrane patch facing this submembrane voxel:
$iCa=iCaL+iCaT+ibCa−2∙iNCX.$
ICaL is included in this equation only for submembrane voxels facing a CRU (ICaL is injected and equally distributed among the nine subspace voxels of the respective dyadic space). The local Ca currents iCaL, iCaT, and ibCa have inward direction and (by convention) are defined as negative. Therefore, the minus sign before iCa in Eq. 5 ensures positive change in [Ca] in the submembrane voxel by respective Ca influx. During diastole local iNCX also flows inwardly, but NCX exchanges 1 Ca ion to 3 Na ions. Hence, the inward iNCX generates a Ca efflux. That is why it has a different sign.
In a voxel of ring layer:
$∂Cacyt,ring∂t=Dcyt∇2Cacyt,ring−jup,ringvFSR,ringvcyt,ring+∑jcyt_sub_ringvcyt,subvcyt,ring−jcyt_ring_to_core−jCM.$
(6)
$∑jcyt_sub_ring$ is the sum of diffusion fluxes from neighboring smaller submembrane voxels
$∂CaFSR,ring∂t=DFSR∇2CaFSR,ring+jup,ring−jFSR_ring_to_JSR−jFSR_ring_to_core.$
(7)

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

In a given JSR:
$∂CaJSR∂t=∑jFSR_ring_to_JSRvring,cytvJSR−Ispark/(2∙F∙vJSR)−jCQ.$
(8)
where Ispark is given by Eq. 3.
jFSR_ring_to_JSR is the sum of diffusion fluxes between JSR and respective FSR parts of the neighboring ring voxels. jCQ is Ca flux of Ca buffering by calsequestrin:
$jCQ=CQtot·∂fCQ∂t.$
(9)
$dfCQ/dt=kfCQ∙CajSR∙(1−fCQ)−kbCQ∙fCQ(t).$
(10)
In the core:
$∂Cacyt,core∂t=∑Jcyt_ring_to_corevcyt,ringvcyt,core−jup,corevFSR,corevcyt,core−jCM.$
(11)
$∂CaFSR,core∂t=∑JFSR_ring_to_corevFSR,ringvFSR,core+jup,core.$
(12)
where ∑jcyt_ring_to_core and ∑jFSR_ring_to_core are the respective sums of diffusion fluxes in cytosol and FSR of all ring voxels.
In any cytoplasmic voxel (subspace, ring, and the core):jCM is Ca flux of Ca buffering by calmodulin:
$jCM=CMtot·∂fCM/∂t.$
(13)
$∂fCM∂t=kfCM∙Cacyt∙(1−fCMs)−kbCM∙fCM(x,y,t).$
(14)

4. Electrophysiology

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 Ist and background Na current IbNa 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 INCX density and Ca-dependent ICaL 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.

4.1 Fixed ion concentrations, mM

Cao = 2: Extracellular Ca concentration; Ko = 5.4: Extracellular K concentration; Ki= 140: Intracellular K concentration; Nao = 140: Extracellular Na concentration; Nai = 10: Intracellular Na concentration.

4.2 The Nernst equation and electric potentials, mV

ENa = ET ∙ ln(Nao/Nai): Equilibrium potential for Na; EK = ET ∙ ln(Ko/Ki): Equilibrium potential for K; EKs= ET∙ ln[(Ko+ 0.12∙ Nao)/(Ki+ 0.12 ∙ Nai)]: Reversal potential of IKs. Where ET is “RT/F” factor = 26.72655 mV at 37°C, ECaL= 45: Apparent reversal potential of ICaL; ECaT= 45: Apparent reversal potential of ICaT.

4.3 Membrane potential, Vm

Net membrane current determines time derivative of the membrane potential.
$dVm/dt=−ICaL+ICaT+IKr+IKs+Ito+Isus+If+INaK+IbCa+INCX/Cm.$

4.4 Formulation of cell membrane ion currents

Kinetics of ion currents are described by gating variables (described below) in respective differential equations:
$dyi/dt=(yi,∞−y)/τyi,$
$(yi=dL,fL,fCa,dT,fT,paF,paS,pi_,n,q,r,y).$
τyi: time constant for a gating variable yi; αyi and βyi: opening and closing rates for channel gating; yi,: steady-state curve for a gating variable yi.

L-type Ca current (ICaL)

The whole-cell ICaL is calculated as a sum of local currents iCaL,i in each CRU.
$ICaL=∑i∈CRUiCaL,i.$
This reflects reports that LCCs are colocalized with RyRs (Christel et al., 2012). Thus, the whole-cell maximum ICaL conductance (gCaL) is distributed locally and equally among all CRUs:
$gCaL,i=gCaL/NCRU.$
The respective local iCaL,i in each CRU is calculated based on formulations of Kurata et al. (2002) and subsequent modifications in our common pool models (Maltsev and Lakatta, 2009; Maltsev and Lakatta, 2010; Maltsev and Lakatta, 2013), but now its Ca-dependent inactivation (fCa,i) is determined by local subspace [Ca] (Casub,i):
$iCaL,i=Cm·gCaL,i·(Vm−ECaL)·dL·fL·fCa,i,$
$dL,∞=1/1+exp−Vm−VdL6,$
$fL,∞=1/{1+exp[(Vm+35)7.3]},$
$αdL=−0.02839·Vm+35/exp−Vm+35/2.5−1−0.0849·Vm/exp−Vm/4.8−1,$
$βdL=0.01143·Vm−5expVm−52.5−1,$
$τdL=1/(αdL+βdL),$
$τfL=k_tau_fL·(257.1·exp{−[Vm+32.513.9]2}+44.3),$
$fCa,∞=KmfCa/(KmfCa+Casub,i),$
$τfCa=fCa,∞/αfCa.$

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

We set the midpoint of ICaL activation VdL to −6.6 mV as in SAN cell models of Wilders et al. (1991) and Dokos et al. (1996). This value of VdL is relatively high with respect to other SAN cell models and therefore was considered in our study to simulate ICaL generated by Cav1.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 Cav1.3 into ICaL. Both isoforms Cav1.2 and Cav1.3 are colocalized with RyRs (Christel et al., 2012), but Cav1.3 has a lower (more hyperpolarized) voltage activation threshold. Therefore, we simulated ICaL generated by Cav1.3 using the same formulations above, but VdL was set to −13.5 mV. The resultant VdL shift increased dL,∞ of Cav1.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 Cav1.3 in gCaL from 0 to 100%, with total conductance gCaL of Cav1.2 and Cav1.3 remaining constant (gCav1.2 + gCav1.3 = gCaL = 0.464 nS/pF = const).

T-type Ca current (ICaT)

It is based on formulations of Demir et al. (1994) and modified by Kurata et al. (2002).
$ICaT=Cm∙gCaT,max∙(Vm−ECaT)∙dT∙fT,$
$dT,∞=1/{1+exp[−(Vm+26.3)/6.0]},$
$fT,∞=1/{1+exp[(Vm+61.7)/5.6]},$
$τdT=1/1.068∙expVm+26.3/30+1.068∙exp−Vm+26.3/30,$
$τfT=1{0.0153∙exp[−(Vm+61.7)83.3]+0.015∙exp[(Vm+61.7)15.38]}.$

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

Rapidly activating delayed rectifier K+current (IKr)

It is based on formulations of Zhang et al. (2000), further modified by Kurata et al. (2002).
$IKr=Cm∙gKr,max∙(Vm−EK)∙(0.6∙paF+0.4∙paS)∙pi,$
$pa,∞=1{1+exp[−(Vm+23.2)10.6]},$
$pi,∞=1{1+exp[(Vm+28.6)17.1]},$
$τpaF=k_tau_IKr∙0.84655354/0.0372∙expVm/15.9+0.00096∙exp−Vm/22.5,$
$τpaS=k_tau_IKr∙0.84655354/[0.0042∙exp(Vm/17.0)+0.00015∙exp(−Vm/21.6)],$
$τpi=1[0.1∙exp(−Vm54.645)+0.656∙exp(Vm106.157)].$

Slowly activating delayed rectifier K+current (IKs)

It is based on formulations of Zhang et al. (2000).
$IKs=Cm∙gKs,max∙(Vm−EKs)∙n2,$
$αn=0.014{1+exp[−(Vm−40)9]},$
$βn=0.001∙exp(−Vm45),$
$n∞=αnαn+βn,$
$τn=1αn+βn.$

4-aminopyridine-sensitive currents (I4AP= Ito+ Isus)

It is based on formulations of Zhang et al. (2000).
$Ito=Cm∙gto,max∙(Vm−EK)∙q∙r,$
$Isus=Cm∙gsus,max∙(Vm−EK)∙r,$
$q∞=1/{1+exp[(Vm+49)/13]},$
$r∞=1{1+exp[−(Vm−19.3)15]},$
$τq=39.102/{0.57∙exp[−0.08∙(Vm+44)]+0.065∙exp[0.1∙(Vm+45.93)]}+6.06,$
$τr=14.40516/1.037∙exp0.09∙Vm+30.61+0.369∙exp−0.12∙Vm+23.84+2.75352.$

Hyperpolarization-activated, “funny” current (If)

It is based on formulations of Wilders et al. (1991) and Kurata et al. (2002).
$If=IfNa+IfK,$
$y∞=1{1+exp[(Vm−VIf,12)13.5]},$
$τy=0.7166529{exp[−(Vm+386.9)45.302]+exp[(Vm−73.08)19.231]},$
$IfNa=Cm∙0.3833∙gIf,max∙(Vm−ENa)∙y2,$
$IfK=Cm∙0.6167∙gIf,max∙(Vm−EK)∙y2.$

Na+-K+pump current (INaK)

It is based on formulations of Kurata et al. (2002), which were in turn based on the experimental work of Sakai et al. (1996) for rabbit SAN cell.
$INaK=Cm∙INaK,max∙1+(KmKp/Ko)1.2−1∙1+(KmNap/Nai)1.3−1∙1+exp−Vm−ENa+120/30−1.$
Ca-background current (IbCa)
$IbCa=Cm∙gbCa∙(Vm−ECaL).$

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

INCX

It is based on original formulations from Dokos et al. (1996). INCX is modulated by local Ca and therefore the whole-cell INCX was calculated as a sum of local currents INCX,i in respective membrane patches facing each submembrane voxel. Thus,
$INCX=∑i=1NvoxelsINCX,i.$
For membrane voltage Vm in each i-th membrane patch with subspace Casub,i, the respective local NCX current (INCX,i) was calculated as follows:
$INCX,i=Cm∙kNCX∙k21∙x2−k12∙x1x1+x2+x3+x4,$
$do=1+CaoKco∙1+expQco∙VmET+NaoK1no∙1+NaoK2no∙1+NaoK3no,$
$k43=NaiK3ni+Nai,$
$k41=exp[−Qn∙Vm/(2ET)],$
$k34=NaoK3no+Nao,$
$k21=(CaoKCo)∙exp(Qco∙VmET)do,$
$k23=(NaoK1no)∙(NaoK2no)∙(1+NaoK3no)∙exp[−Qn∙Vm2ET]do,$
$k32=exp[Qn∙Vm/(2ET)],$
$x1=k34∙k41∙k23+k21+k21∙k32∙k43+k41,$
$di=1+Casub,iKci∙1+exp−Qci∙VmET+NaiKcni+NaiK1ni∙1+NaiK2ni∙1+NaiK3ni,$
$k12=(Casub,i/Kci)∙exp(−Qci∙Vm/ET)/di,$
$k14=(Nai/K1ni)∙(Nai/K2ni)∙(1+Nai/K3ni)∙exp[Qn∙Vm/(2ET)]/di,$
$x2=k43∙k32∙(k14+k12)+k41∙k12∙(k34+k32),$
$x3=k43∙k14∙(k23+k21)+k12∙k23∙(k43+k41),$
$x4=k34∙k23∙(k14+k12)+k21∙k14∙(k34+k32).$

5. Initial values

Initial values for electrophysiology are in Table A1.

Initial values for Ca dynamics are in Table A2.

6. Summary of model parameters

6.1 Fixed ion concentrations

Cao = 2 mM: Extracellular [Ca]; Ko = 5.4 mM: Extracellular [K]; Nao = 140 mM: Extracellular [Na]; Ki = 140 mM: Intracellular [K]; Nai = 10 mM: Intracellular [Na].

6.2 Membrane currents

ECaL = 45 mV: Apparent reversal potential of ICaL; gCaL = 0.464 nS/pF: Conductance of ICaL; KmfCa = 0.03 mM: Dissociation constant of Ca-dependent ICaL inactivation; betafCa = 60 mM−1 · ms−1: Ca association rate constant for ICaL; alfafCa = 0.021 ms−1: Ca dissociation rate constant for ICaL; k_tau_fL = 0.5: Scaling factor for tau_fL used to tune the model; ECaT = 45 mV: Apparent reversal potential of ICaT, mV; gCaT = 0.1832 nS/pF: Conductance of ICaT; gIf = 0.105 nS/pF: Conductance of If; VIf,1/2 = −64 mV: half activation voltage of If; gKr = 0.05679781 nS/pF: Conductance of delayed rectifier K current rapid component; k_tau_IKr = 0.3: scaling factor for tau_paF and tau_paS used to tune the model; gKs = 0.0259 nS/pF: Conductance of delayed rectifier K current slow component; gto = 0.252 nS/pF: Conductance of 4-aminopyridine sensitive transient K current; gsus = 0.02 nS/pF: Conductance of 4-aminopyridine sensitive sustained K current; INaKmax = 1.44 pA/pF: Maximum Na+/K+ pump current; KmKp = 1.4 mM: Half-maximal Ko for INaK; KmNap = 14 mM: Half-maximal Nai for INaK; gbCa = 0.003 nS/pF: Conductance of background Ca current; kNCX = 48.75 pA/pF: Maximum amplitude of INCX.

6.3 Dissociation constants for NCX

K1ni = 395.3 mM: intracellular Na binding to first site on NCX; K2ni = 2.289 mM: intracellular Na binding to second site on NCX; K3ni = 26.44 mM: intracellular Na binding to third site on NCX; K1no = 1628 mM: extracellular Na binding to first site on NCX; K2no = 561.4 mM: extracellular Na binding to second site on NCX; K3no = 4.663 mM: extracellular Na binding to third site on NCX; Kci = 0.0207 mM: intracellular Ca binding to NCX transporter; Kco = 3.663 mM: extracellular Ca binding to NCX transporter; Kcni = 26.44 mM: intracellular Na and Ca simultaneous binding to NCX.

6.4 NCX fractional charge movement

Qci = 0.1369: intracellular Ca occlusion reaction of NCX; Qco = 0: extracellular Ca occlusion reaction of NCX; Qn = 0.4315: Na occlusion reactions of NCX.

6.5 Ca buffering

kbCM = 0.542 ms−1: Ca dissociation constant for calmodulin; kfCM = 227.7 mM−1· ms−1: Ca association constant for calmodulin; kbCQ = 0.445 ms−1: Ca dissociation constant for calsequestrin; kfCQ = 0.534 mM−1· ms−1: Ca association constant for calsequestrin; CQtot = 30 mM: Total calsequestrin concentration; CMtot = 0.045 mM: Total calmodulin concentration.

6.6 SR Ca ATPase function

Kmf = 0.000246 mM: the cytosolic side Kd of SR Ca pump; Kmr = 1.7 mM: the lumenal side Kd of SR Ca pump; H = 1.787: cooperativity of SR Ca pump; Pup = 0.014 mM/ms: Maximal rate of Ca uptake by SR Ca pump.

6.7 CRU (Ca release and JSR)

CRU_Casens = 0.00015 mM: sensitivity of Ca release to Casub; CRU_ProbConst = 0.00027 ms−1: CRU open probability rate at Casub = CRU_Casens; CRU_ProbPower = 3: Cooperativity of CRU activation by Casub; 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.

6.8 Cell geometry, compartments, and voxels

Lcell = 53.28 μm: Cell length; rcell = 3.437747 μm: Cell radius; Cm = 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.

6.9 Ca diffusion

Dcyt = 0.35 μm2/ms: Diffusion coefficient of free Ca in cytosol; DFSR = 0.06 μm2/ms: Diffusion coefficient of free Ca in FSR.

7. gCaL sensitivity analysis

For square lattice and uniform random distributions of CRUs, we performed sensitivity analysis for gCaL 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 Vm time series are given Fig. S2 for each CRU distribution.

8. Simulations of βAR stimulation effect

Effect of βAR stimulation was modeled essentially as we previously reported (Maltsev and Lakatta, 2010) by increasing ICaL, IKr, If, and Ca uptake rate by FSR via SERCA pumping. Specifically, whole-cell maximum ICaL conductance gCaL was increased by a factor of 1.75 from 0.464 to 0.812 nS/pF; whole-cell maximum IKr conductance gKr was increased by a factor of 1.5 from 0.05679781 to 0.085196715 nS/pF; the midpoint of If activation curve was shifted to more depolarized potential by 7.8 mV from −64 to −56.2 mV; and the maximum Ca uptake rate Pup 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.

9. Model integration

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 Cav1.2 and Cav1.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.

Asghari
,
P.
,
D.R.
Scriven
,
M.
Ng
,
P.
Panwar
,
K.C.
Chou
,
F.
van Petegem
, and
E.D.
Moore
.
2020
.
Cardiac ryanodine receptor distribution is dynamic and changed by auxiliary proteins and post-translational modification
.
Elife
.
9
:e51602.
Baig
,
S.M.
,
A.
Koschak
,
A.
Lieb
,
M.
Gebhart
,
C.
Dafinger
,
G.
Nurnberg
,
A.
Ali
,
I.
,
M.J.
Sinnegger-Brauns
,
N.
Brandt
, et al
.
2011
.
Loss of Ca(v)1.3 (CACNA1D) function in a human channelopathy with bradycardia and congenital deafness
.
Nat. Neurosci.
14
:
77
84
.
Bogdanov
,
K.Y.
,
T.M.
, and
E.G.
Lakatta
.
2001
.
Sinoatrial nodal cell ryanodine receptor and Na+-Ca2+ exchanger: Molecular partners in pacemaker regulation
.
Circ. Res.
88
:
1254
1258
.
Brennan
,
J.A.
,
Q.
Chen
,
A.
Gams
,
J.
Dyavanapalli
,
D.
Mendelowitz
,
W.
Peng
, and
I.R.
Efimov
.
2020
.
Evidence of superior and inferior sinoatrial nodes in the mammalian heart
.
JACC. Clin. Electrophysiol.
6
:
1827
1840
.
Bychkov
,
R.
,
M.
Juhaszova
,
K.
Tsutsui
,
C.
Coletta
,
M.D.
Stern
,
V.A.
Maltsev
, and
E.G.
Lakatta
.
2020
.
Synchronized cardiac impulses emerge from heterogeneous local calcium signals within and among cells of pacemaker tissue
.
JACC Clin. Electrophysiol.
6
:
907
931
.
Chen-Izu
,
Y.
,
C.W.
Ward
,
W.
Stark
Jr.
,
T.
Banyasz
,
M.P.
Sumandea
,
C.W.
Balke
,
L.T.
Izu
, and
X.H.
Wehrens
.
2007
.
Phosphorylation of RyR2 and shortening of RyR2 cluster spacing in spontaneously hypertensive rat with heart failure
.
Am. J. Physiol. Heart Circ. Physiol.
293
:
H2409
H2417
.
Chen
,
B.
,
Y.
Wu
,
P.J.
Mohler
,
M.E.
Anderson
, and
L.S.
Song
.
2009
.
Local control of Ca2+-induced Ca2+ release in mouse sinoatrial node cells
.
J. Mol. Cell Cardiol.
47
:
706
715
.
Cheng
,
H.
,
W.J.
Lederer
, and
M.B.
Cannell
.
1993
.
Calcium sparks: Elementary events underlying excitation-contraction coupling in heart muscle
.
Science
.
262
:
740
744
.
Christel
,
C.J.
,
N.
Cardona
,
P.
Mesirca
,
S.
Herrmann
,
F.
Hofmann
,
J.
Striessnig
,
A.
Ludwig
,
M.E.
Mangoni
, and
A.
Lee
.
2012
.
Distinct localization and modulation of Cav1.2 and Cav1.3 L-type Ca2+ channels in mouse sinoatrial node
.
J. Physiol
.
590
:
6327
6342
.
Clancy
,
C.E.
, and
L.F.
Santana
.
2020
.
Evolving discovery of the origin of the heartbeat: A new perspective on sinus rhythm
.
JACC. Clin. Electrophysiol.
6
:
932
934
.
Demir
,
S.S.
,
J.W.
Clark
,
C.R.
Murphey
, and
W.R.
Giles
.
1994
.
A mathematical model of a rabbit sinoatrial node cell
.
Am. J. Physiol.
266
:
C832
C852
.
DiFrancesco
,
D.
, and
C.
Tromba
.
1988
.
Inhibition of the hyperpolarization-activated current (if) induced by acetylcholine in rabbit sino-atrial node myocytes
.
J. Physiol.
405
:
477
491
.
Dokos
,
S.
,
B.
Celler
, and
N.
Lovell
.
1996
.
Ion currents underlying sinoatrial node pacemaker activity: A new single cell mathematical model
.
J. Theor. Biol.
181
:
245
272
.
Fenske
,
S.
,
K.
Hennis
,
R.D.
Rotzer
,
V.F.
Brox
,
E.
Becirovic
,
A.
Scharr
,
C.
Gruner
,
T.
Ziegler
,
V.
Mehlfeld
,
J.
Brennan
, et al
.
2020
.
cAMP-dependent regulation of HCN4 controls the tonic entrainment process in sinoatrial node pacemaker cells
.
Nat. Commun.
11
:
5555
.
Gratz
,
D.
,
B.
Onal
,
A.
Dalic
, and
T.J.
Hund
.
2018
.
Synchronization of pacemaking in the sinoatrial node: A mathematical modeling study
.
Front. Phys.
6
:
63
.
Greiser
,
M.
,
H.C.
Joca
, and
W.J.
Lederer
.
2019
.
Pacemaker organization at the nanoscale: Imaging of ryanodine receptors as clusters in single sinoatrial nodal cells
.
Biophys. J.
116
:
380A
.
Guarina
,
L.
,
A.N.
Moghbel
,
M.S.
,
R.H.
Cudmore
,
D.
Sato
,
C.E.
Clancy
, and
L.F.
Santana
.
2022
.
Biological noise is a key determinant of the reproducibility and adaptability of cardiac pacemaking and EC coupling
.
J. Gen. Physiol.
154
:e202012613.
Himeno
,
Y.
,
N.
Sarai
,
S.
Matsuoka
, and
A.
Noma
.
2008
.
Ionic mechanisms underlying the positive chronotropy induced by beta1-adrenergic stimulation in Guinea pig sinoatrial node cells: A simulation study
.
J. Physiol. Sci.
58
:
53
65
.
Honjo
,
H.
,
M.R.
Boyett
,
I.
Kodama
, and
J.
Toyama
.
1996
.
Correlation between electrical activity and the size of rabbit sino-atrial node cells
.
J. Physiol.
496 ( Pt 3)
:
795
808
.
Huser
,
J.
,
L.A.
Blatter
, and
S.L.
Lipsius
.
2000
.
Intracellular Ca2+ release contributes to automaticity in cat atrial pacemaker cells
.
J. Physiol.
524 Pt 2
:
415
422
.
Imtiaz
,
M.S.
,
P.Y.
von der Weid
,
D.R.
Laver
, and
D.F.
van Helden
.
2010
.
SR Ca2+ store refill: A key factor in cardiac pacemaking
.
J. Mol. Cell. Cardiol.
49
:
412
426
.
,
S.
,
H.
Zhang
,
J.O.
Tellez
,
N.
Shibata
,
K.
Nakazawa
,
K.
Kamiya
,
I.
Kodama
,
K.
Mitsui
,
H.
Dobrzynski
,
M.R.
Boyett
, and
H.
Honjo
.
2014
.
Importance of gradients in membrane properties and electrical coupling in sinoatrial node pacing
.
PLoS One
.
9
:e94565.
Kim
,
M.S.
,
A.V.
Maltsev
,
O.
Monfredi
,
L.A.
Maltseva
,
A.
Wirth
,
M.C.
Florio
,
K.
Tsutsui
,
D.R.
Riordon
,
S.P.
Parsons
,
S.
Tagirova
, et al
.
2018
.
Heterogeneity of calcium clock functions in dormant, dysrhythmically and rhythmically firing single pacemaker cells isolated from SA node
.
Cell Calcium
.
74
:
168
179
.
Kim
,
M.S.
,
O.
Monfredi
,
L.A.
Maltseva
,
E.G.
Lakatta
, and
V.A.
Maltsev
.
2021
.
Beta-adrenergic stimulation synchronizes a broad spectrum of action potential firing rates of cardiac pacemaker cells toward a higher population average
.
Cells
.
10
:
2124
.
Kurata
,
Y.
,
I.
Hisatome
,
S.
Imanishi
, and
T.
Shibamoto
.
2002
.
Dynamical description of sinoatrial node pacemaking: Improved mathematical model for primary pacemaker cell
.
Am. J. Physiol. Heart Circ. Physiol.
283
:
H2074
H2101
.
Kurata
,
Y.
,
I.
Hisatome
, and
T.
Shibamoto
.
2012
.
Roles of sarcoplasmic reticulum Ca2+ cycling and Na+/Ca2+ exchanger in sinoatrial node pacemaking: Insights from bifurcation analysis of mathematical models
.
Am. J. Physiol. Heart Circ. Physiol.
302
:
H2285
H2300
.
Lakatta
,
E.G.
,
V.A.
Maltsev
, and
T.M.
.
2010
.
A coupled SYSTEM of intracellular Ca2+ clocks and surface membrane voltage clocks controls the timekeeping mechanism of the heart’s pacemaker
.
Circ. Res.
106
:
659
673
.
Laver
,
D.R.
,
C.H.
Kong
,
M.S.
Imtiaz
, and
M.B.
Cannell
.
2013
.
Termination of calcium-induced calcium release by induction decay: An emergent property of stochastic channel gating and molecular scale architecture
.
J. Mol. Cell Cardiol.
54
:
98
100
.
Li
,
K.
,
Z.
Chu
, and
X.
Huang
.
2018
.
Annihilation of the pacemaking activity in the sinoatrial node cell and tissue
.
8
:
125319
.
,
J.
,
O.
Bortolotti
,
E.
Torre
,
I.
Bidaud
,
N.
Lamb
,
A.
Fernandez
,
J.Y.
Le Guennec
,
M.E.
Mangoni
, and
P.
Mesirca
.
2022
.
L-type Cav1.3 calcium channels are required for beta-adrenergic triggered automaticity in dormant mouse sinoatrial pacemaker cells
.
Cells
.
11
:
1114
.
Lyashkov
,
A.E.
,
J.
Behar
,
E.G.
Lakatta
,
Y.
Yaniv
, and
V.A.
Maltsev
.
2018
.
Positive feedback mechanisms among local Ca releases, NCX, and ICaL ignite pacemaker action potentials
.
Biophys. J.
114
:
2024
.
Lyashkov
,
A.E.
,
M.
Juhaszova
,
H.
Dobrzynski
,
T.M.
,
V.A.
Maltsev
,
O.
Juhasz
,
H.A.
Spurgeon
,
S.J.
Sollott
, and
E.G.
Lakatta
.
2007
.
Calcium cycling protein density and functional importance to automaticity of isolated sinoatrial nodal cells are independent of cell size
.
Circ. Res.
100
:
1723
1731
.
Lyashkov
,
A.E.
,
T.M.
,
I.
Zahanich
,
Y.
Li
,
A.
Younes
,
H.B.
Nuss
,
H.A.
Spurgeon
,
V.A.
Maltsev
, and
E.G.
Lakatta
.
2009
.
Cholinergic receptor signaling modulates spontaneous firing of sinoatrial nodal cells via integrated effects on PKA-dependent Ca2+ cycling and IKACh
.
Am. J. Physiol. Heart Circ. Physiol.
297
:
H949
H959
.
Maltsev
,
A.V.
,
V.A.
Maltsev
,
M.
Mikheev
,
L.A.
Maltseva
,
S.G.
Sirenko
,
E.G.
Lakatta
, and
M.D.
Stern
.
2011
.
Synchronization of stochastic Ca2+ release units creates a rhythmic Ca2+ clock in cardiac pacemaker cells
.
Biophys. J.
100
:
271
283
.
Maltsev
,
A.V.
,
V.A.
Maltsev
, and
M.D.
Stern
.
2017a
.
Clusters of calcium release channels harness the Ising phase transition to confine their elementary intracellular signals
.
.
114
:
7525
7530
.
Maltsev
,
A.V.
,
V.A.
Maltsev
, and
M.D.
Stern
.
2017b
.
Stabilization of diastolic calcium signal via calcium pump regulation of complex local calcium releases and transient decay in a computational model of cardiac pacemaker cell with individual release channels
.
PLoS Comput. Biol.
13
:e1005675.
Maltsev
,
A.V.
,
M.D.
Stern
,
E.G.
Lakatta
, and
V.A.
Maltsev
.
2022
.
Functional heterogeneity of cell populations increases robustness of pacemaker function in a numerical model of the sinoatrial node tissue
.
Front. Physiol.
13
:
845634
.
Maltsev
,
A.V.
,
M.D.
Stern
, and
V.A.
Maltsev
.
2019
.
Mechanisms of calcium leak from cardiac sarcoplasmic reticulum revealed by statistical mechanics
.
Biophys. J.
116
:
2212
2223
.
Maltsev
,
A.V.
,
Y.
Yaniv
,
M.D.
Stern
,
E.G.
Lakatta
, and
V.A.
Maltsev
.
2013
.
RyR-NCX-SERCA local crosstalk ensures pacemaker cell function at rest and during the fight-or-flight reflex
.
Circ. Res.
113
:
e94
e100
.
Maltsev
,
V.A.
, and
E.G.
Lakatta
.
2009
.
Synergism of coupled subsarcolemmal Ca2+ clocks and sarcolemmal voltage clocks confers robust and flexible pacemaker function in a novel pacemaker cell model
.
Am. J. Physiol. Heart Circ. Physiol.
296
:
H594
H615
.
Maltsev
,
V.A.
, and
E.G.
Lakatta
.
2010
.
A novel quantitative explanation for the autonomic modulation of cardiac pacemaker cell automaticity via a dynamic system of sarcolemmal and intracellular proteins
.
Am. J. Physiol. Heart Circ. Physiol.
298
:
H2010
H2023
.
Maltsev
,
V.A.
, and
E.G.
Lakatta
.
2013
.
Numerical models based on a minimal set of sarcolemmal electrogenic proteins and an intracellular Ca(2+) clock generate robust, flexible, and energy-efficient cardiac pacemaking
.
J. Mol. Cell Cardiol.
59
:
181
195
.
Maltsev
,
V.A.
,
A.V.
Maltsev
,
M.
Juhaszova
,
S.
Sirenko
,
O.
Monfredi
,
H.
Shroff
,
A.
York
,
S.J.
Sollott
,
E.G.
Lakatta
, and
M.D.
Stern
.
2016
.
Cardiac pacemaker cell function at a super-resolution scale of SIM: Distribution of RyRs, calcium dynamics, and numerical modeling
.
Biophys. J.
110
:
267A
.
Maltsev
,
V.A.
,
Y.
Yaniv
,
A.V.
Maltsev
,
M.D.
Stern
, and
E.G.
Lakatta
.
2014
.
Modern perspectives on numerical modeling of cardiac pacemaker cell
.
J. Pharmacol. Sci.
125
:
6
38
.
McDonnell
,
M.D.
, and
D.
Abbott
.
2009
.
What is stochastic resonance? Definitions, misconceptions, debates, and its relevance to biology
.
PLoS Comput. Biol.
5
:e1000348
Mesirca
,
P.
,
A.G.
Torrente
, and
M.E.
Mangoni
.
2015
.
Functional role of voltage gated Ca2+ channels in heart automaticity
.
Front. Physiol.
6
:
19
.
Monfredi
,
O.
,
K.
Tsutsui
,
B.
Ziman
,
M.D.
Stern
,
E.G.
Lakatta
, and
V.A.
Maltsev
.
2018
.
Electrophysiological heterogeneity of pacemaker cells in rabbit intercaval region, including SA node: Insights from recording multiple ion currents in each cell
.
Am. J. Physiol. Heart Circ. Physiol.
314
:
H403
H414
.
Musa
,
H.
,
M.
Lei
,
H.
Honjo
,
S.A.
Jones
,
H.
Dobrzynski
,
M.K.
Lancaster
,
Y.
Takagishi
,
Z.
Henderson
,
I.
Kodama
, and
M.R.
Boyett
.
2002
.
Heterogeneous expression of Ca2+ handling proteins in rabbit sinoatrial node
.
J. Histochem. Cytochem.
50
:
311
324
.
Nivala
,
M.
,
C.Y.
Ko
,
M.
Nivala
,
J.N.
Weiss
, and
Z.
Qu
.
2012
.
Criticality in intracellular calcium signaling in cardiac myocytes
.
Biophys. J.
102
:
2433
2442
.
Oren
,
R.V.
, and
C.E.
Clancy
.
2010
.
Determinants of heterogeneity, excitation and conduction in the sinoatrial node: A model study
.
Plos Comput. Biol.
6
:e1001041.
Qu
,
Z.
,
A.
Garfinkel
,
J.N.
Weiss
, and
M.
Nivala
.
2011
.
Multi-scale modeling in biology: How to bridge the gaps between scales?
.
Prog. Biophys. Mol. Biol.
107
:
21
31
.
Rigg
,
L.
,
B.M.
Heath
,
Y.
Cui
, and
D.A.
Terrar
.
2000
.
Localisation and functional significance of ryanodine receptors during beta-adrenoceptor stimulation in the guinea-pig sino-atrial node
.
Cardiovasc. Res.
48
:
254
264
.
Sakai
,
R.
,
N.
Hagiwara
,
N.
Matsuda
,
H.
Kassanuki
, and
S.
Hosoda
.
1996
.
Sodium--potassium pump current in rabbit sino-atrial node cells
.
J. Physiol.
490 (Pt 1)
:
51
62
.
Severi
,
S.
,
M.
Fantini
,
L.A.
Charawi
, and
D.
DiFrancesco
.
2012
.
An updated computational model of rabbit sinoatrial action potential to investigate the mechanisms of heart rate modulation
.
J. Physiol.
590
:
4483
4499
.
Shannon
,
T.R.
,
F.
Wang
,
J.
Puglisi
,
C.
Weber
, and
D.M.
Bers
.
2004
.
A mathematical treatment of integrated Ca dynamics within the ventricular myocyte
.
Biophys. J.
87
:
3351
3371
.
Stern
,
M.D.
1992
.
Theory of excitation-contraction coupling in cardiac muscle
.
Biophys. J.
63
:
497
517
.
Stern
,
M.D.
,
L.A.
Maltseva
,
M.
Juhaszova
,
S.J.
Sollott
,
E.G.
Lakatta
, and
V.A.
Maltsev
.
2014
.
Hierarchical clustering of ryanodine receptors enables emergence of a calcium clock in sinoatrial node cells
.
J. Gen. Physiol.
143
:
577
604
.
Stern
,
M.D.
,
G.
Pizarro
, and
E.
Rios
.
1997
.
Local control model of excitation-contraction coupling in skeletal muscle
.
J. Gen. Physiol.
110
:
415
440
.
Stern
,
M.D.
,
E.
Rios
, and
V.A.
Maltsev
.
2013
.
Life and death of a cardiac calcium spark
.
J. Gen. Physiol.
142
:
257
274
.
Ter Keurs
,
H.E.D.J.
, and
P.A.
Boyden
.
2007
.
Calcium and arrhythmogenesis
.
Physiol. Rev.
87
:
457
506
.
Torrente
,
A.G.
,
P.
Mesirca
,
P.
Neco
,
R.
Rizzetto
,
S.
Dubel
,
C.
Barrere
,
M.
Sinegger-Brauns
,
J.
Striessnig
,
S.
Richard
,
J.
Nargeot
, et al
.
2016
.
L-type Cav1.3 channels regulate ryanodine receptor-dependent Ca2+ release during sino-atrial node pacemaker activity
.
Cardiovasc. Res
.
109
:
451
461
.
Tsutsui
,
K.
,
M.C.
Florio
,
A.
Yang
,
A.N.
Wirth
,
D.
Yang
,
M.S.
Kim
,
B.D.
Ziman
,
R.
Bychkov
,
O.J.
Monfredi
,
V.A.
Maltsev
, and
E.G.
Lakatta
.
2021
.
cAMP-dependent signaling restores AP firing in dormant SA node cells via enhancement of surface membrane currents and calcium coupling
.
Front. Physiol.
12
:
596832
.
Tsutsui
,
K.
,
O.
Monfredi
, and
E.G.
Lakatta
.
2016
.
A general theory to explain heart rate and cardiac contractility changes with age
.
Physiol. Mini Rev.
9
:
9
25
Tsutsui
,
K.
,
O.J.
Monfredi
,
S.G.
Sirenko-Tagirova
,
L.A.
Maltseva
,
R.
Bychkov
,
M.S.
Kim
,
B.D.
Ziman
,
K.V.
Tarasov
,
Y.S.
Tarasova
,
J.
Zhang
, et al
.
2018
.
A coupled-clock system drives the automaticity of human sinoatrial nodal pacemaker cells
.
Sci. Signal.
11
:eaap7608.
Veron
,
G.
,
V.A.
Maltsev
,
M.D.
Stern
, and
A.V.
Maltsev
.
2021
.
Elementary intracellular signals are initiated by a transition of release channel system from a metastable state
.
Arxiv
.
,
T.M.
,
K.Y.
Bogdanov
, and
E.G.
Lakatta
.
2002
.
beta-Adrenergic stimulation modulates ryanodine receptor Ca2+ release during diastolic depolarization to accelerate pacemaker activity in rabbit sinoatrial nodal cells
.
Circ. Res.
90
:
73
79
.
,
T.M.
,
D.X.
Brochet
,
S.
Sirenko
,
Y.
Li
,
H.
Spurgeon
, and
E.G.
Lakatta
.
2010
.
Sarcoplasmic reticulum Ca2+ pumping kinetics regulates timing of local Ca2+ releases and spontaneous beating rate of rabbit sinoatrial node pacemaker cells
.
Circ. Res.
107
:
767
775
.
,
T.M.
,
A.E.
Lyashkov
,
W.
Zhu
,
A.M.
Ruknudin
,
S.
Sirenko
,
D.
Yang
,
S.
Deo
,
M.
Barlow
,
S.
Johnson
,
J.L.
Caffrey
, et al
.
2006
.
High basal protein kinase A-dependent phosphorylation drives rhythmic internal Ca2+ store oscillations and spontaneous beating of cardiac pacemaker cells
.
Circ. Res.
98
:
505
514
.
Weiss
,
J.N.
, and
Z.
Qu
.
2020
.
The sinus node: Still mysterious after all these years
.
JACC Clin. Electrophysiol.
6
:
1841
1843
.
Wilders
,
R.
,
H.J.
Jongsma
, and
A.C.
van Ginneken
.
1991
.
Pacemaker activity of the rabbit sinoatrial node. A comparison of mathematical models
.
Biophys. J.
60
:
1202
1216
.
Yuan
,
X.
,
L.N.
Ratajczyk
,
F.
,
H.H.
Valdivia
,
A.V.
Glukhov
, and
D.
Lang
.
2020
.
Hierarchical pacemaker clustering within the rabbit sinoatrial node is driven by dynamic interaction between the components of the coupled-clock system
.
Biophys. J.
118
:
345A
.
Zhang
,
H.
,
A.V.
Holden
,
I.
Kodama
,
H.
Honjo
,
M.
Lei
,
T.
Varghese
, and
M.R.
Boyett
.
2000
.
Mathematical models of action potentials in the periphery and center of the rabbit sinoatrial node
.
Am. J. Physiol. Heart Circ. Physiol.
279
:
H397
H421
.
Zhang
,
Y.
,
J.L.
Ocampo-Espindola
,
I.Z.
Kiss
, and
A.E.
Motter
.
2021
.
.
.
118
:e2024299118.
Zhou
,
P.
,
Y.T.
Zhao
,
Y.B.
Guo
,
S.M.
Xu
,
S.H.
Bai
,
E.G.
Lakatta
,
H.
Cheng
,
X.M.
Hao
, and
S.Q.
Wang
.
2009
.
Beta-adrenergic signaling accelerates and synchronizes cardiac ryanodine receptor response to a single L-type Ca2+ channel
.
.
106
:
18028
18033
.
Zima
,
A.V.
,
E.
Bovo
,
D.M.
Bers
, and
L.A.
Blatter
.
2010
.
Ca2+ spark-dependent and -independent sarcoplasmic reticulum Ca2+ leak in normal and failing rabbit ventricular myocytes
.
J. Physiol.
588
:
4743
4757
.
Zima
,
A.V.
,
E.
Picht
,
D.M.
Bers
, and
L.A.
Blatter
.
2008
.
Termination of cardiac Ca2+ sparks: Role of intra-SR [Ca2+], release flux, and intra-SR Ca2+ diffusion
.
Circ. Res.
103
:
e105
e115
.

This work is part of a special issue on excitation–contraction coupling.