Ca2+ and Membrane Potential Transitions During Action Potentials Are Self-Similar to Each Other and to Variability of AP Firing Intervals Across the Broad Physiologic Range of AP Intervals During Autonomic Receptor Stimulation

Ca²⁺ and V_m transitions occurring throughout action potential (AP) cycles in sinoatrial nodal (SAN) cells are cues that (1) not only regulate activation states of molecules operating within criticality (Ca²⁺ domain) and limit-cycle (V_m domain) mechanisms of a coupled-clock system that underlies SAN cell automaticity, (2) but are also regulated by the activation states of the clock molecules they regulate. In other terms, these cues are both causes and effects of clock molecular activation (recursion). Recently, we demonstrated that Ca²⁺ and V_m transitions during AP cycles in single SAN cells isolated from mice, guinea pigs, rabbits, and humans are self-similar (obey a power law) and are also self-similar to trans-species AP firing intervals (APFIs) of these cells in vitro, to heart rate in vivo, and to body mass. Neurotransmitter stimulation of β-adrenergic receptor or cholinergic receptor–initiated signaling in SAN cells modulates their AP firing rate and rhythm by impacting on the degree to which SAN clocks couple to each other, creating the broad physiologic range of SAN cell mean APFIs and firing interval variabilities. Here we show that Ca²⁺ and V_m domain kinetic transitions (time to AP ignition in diastole and 90% AP recovery) occurring within given AP, the mean APFIs, and APFI variabilities within the time series of APs in 230 individual SAN cells are self-similar (obey power laws). In other terms, these long-range correlations inform on self-similar distributions of order among SAN cells across the entire broad physiologic range of SAN APFIs, regardless of whether autonomic receptors of these cells are stimulated or not and regardless of the type (adrenergic or cholinergic) of autonomic receptor stimulation. These long-range correlations among distributions of Ca²⁺ and V_m kinetic functions that regulate SAN cell clock coupling during each AP cycle in different individual, isolated SAN cells not in contact with each other. Our numerical model simulations further extended our perspectives to the molecular scale and demonstrated that many ion currents also behave self-similar across autonomic states. Thus, to ensure rapid flexibility of AP firing rates in response to different types and degrees of autonomic input, nature “did not reinvent molecular wheels within the coupled-clock system of pacemaker cells,” but differentially engaged or scaled the kinetics of gears that regulate the rate and rhythm at which the “wheels spin” in a given autonomic input context.

Introduction

The heart is a central player within a hierarchical system of clocks operating within the autonomic neurovisceral axis that creates and synchronizes rhythmic functions ranging from milliseconds to days and beyond (Shivkumar et al., 2016; Lakatta, 2021). The heart’s beating rate and rhythm are regulated by autonomic input to sinoatrial nodal (SAN) pacemaker cells that modulates functions within a coupled-clock system intrinsic to SAN cells (Lakatta et al., 2010).

What Is the Coupled-Clock System Within Pacemaker Cells and How Do Clocks Couple to Each Other?

The SAN cell coupled-clock system comprised a calcium “clock,” the sarcoplasmic reticulum (SR), which continuously oscillates Ca²⁺ via a criticality mechanism (Nivala et al., 2012) and phase-like transitions (Maltsev et al., 2011); the Ca²⁺ clock is continuously but variably coupled to a “membrane clock,” an ensemble of surface membrane ion channels that generates current oscillations via a limit-cycle mechanism (Weiss and Qu, 2020). The criticality mechanisms, in turn, are governed by power law and self-similarity across wide scales (Bak, 1999). The “biochemical engine” of the coupled-clock system is a constitutively active, Ca²⁺ calmodulin-dependent adenylyl cyclase (AC) that generates cyclic AMP (cAMP), leading to modulation of cAMP-gated ion channels, exchange protein directly activated by cAMP (EPAC) signaling, and protein kinase A (PKA) and Calcium/calmodulin-dependent protein kinase II (CaMKII)-dependent kinase activities, mechanisms that regulate intracellular Ca²⁺ levels, Ca²⁺ dynamics and membrane potential within SAN cells (Lakatta et al., 2003, 2006, 2008, 2010; Maltsev and Lakatta, 2008; Yaniv et al., 2015). Variable rates and rhythms at which SAN cells fire action potentials (APs) are controlled by the kinetics of subcellular and cell-wide transitions in [Ca²⁺] gradients and the membrane potential (V_m), and the extent to which V_m and Ca²⁺ become coupled during AP cycles in any given epoch. For more details, see Supplementary Discussion.

The well-known variability of AP firing intervals (APFIs) of isolated SAN cells in vitro, of SAN tissue ex vivo, or of heartbeat intervals in vivo (Monfredi et al., 2014; Yaniv et al., 2014a) indicates that coupled-clock system Ca²⁺ and V_m functions during AP cycles never achieve a true steady state from one AP to the next.

These time-dependent Ca²⁺ and V_m domain transitions during APs are cues, which not only regulate activation states of clock molecules, but are also regulated by the activation status of the very molecules they regulate. In other terms, changes in these cues cause changes in clock molecule activation that feed back to change the characteristics of activation cues. This recursive dynamic imparts robustness to SAN cell automaticity (Maltsev and Lakatta, 2009; Lyashkov et al., 2018). The variability in the degree to which Ca²⁺ and membrane clock molecules couple to each other throughout AP cycles is due to transitions (changes) that occur in Ca²⁺ and V_m domain cues throughout AP cycles (Monfredi et al., 2013; Yaniv et al., 2014b).

Spontaneous transitions in subcellular Ca²⁺ and V_m domains that emerge during the spontaneous diastolic depolarization (DD) phase of an AP cycle have been conceptualized as the AP “ignition phase” (Lyashkov et al., 2018). The ignition process in the Ca²⁺ domain is linked to the emergence of local spontaneous, diastolic oscillatory RyR activation, which generates local Ca²⁺ releases (LCRs) that self-organize to form Ca²⁺ wavelets that propagate locally (Bogdanov et al., 2001; Vinogradova et al., 2004). Ca²⁺-dependent activation of the surface membrane electrogenic Na⁺/Ca²⁺ exchanger generates inward current that accelerates diastolic V_m depolarization and couples the clocks. The time at which the rate of this feed-forward crosstalk acutely accelerates to 0.15 V/s marks the onset of the coupled-clock ignition process (Lyashkov et al., 2018).

Following ignition onset, the extent to which the Ca²⁺ and V_m clock become coupled continues to increase throughout the diastolic period as LCRs and Ca²⁺ wavelets emerge at remote areas across the cell and continue to self-organize in time throughout the cellular space, creating an explosive ensemble Ca²⁺ signal that progressively depolarizes the cell membrane; i.e., clock coupling progressively increases. This Ca²⁺-induced change in V_m increases in clock-coupling cues the activation of low-voltage activated Ca²⁺ channels (Ca_v1.3 and Ca_v3.1), resulting in Ca²⁺ influx that contributes to the further organization of the ensemble LCR Ca²⁺ signal via feed-forward electrochemical (Ca²⁺–V_m–Ca²⁺) signaling, when the diastolic V_m enters a range that cues the activation of L-type Ca²⁺ channels (Ca_v1.2). The ignition phase of the coupled-Ca²⁺ and V_m domain subcellular kinetic transitions culminates in the generation of cell-wide events; a marked transition in the rate of V_m depolarization, due to the activation of Ca_v1.2, results in the rapid AP upstroke and Ca²⁺ influx, which, via Ca²⁺-induced Ca²⁺ release from the SR via RyRs, generates an AP-induced cytosolic Ca²⁺ transient (CaT). In other terms, spontaneous, cell-wide Ca²⁺ signals and APs in SAN cells emerge from spatiotemporal self-organization of spontaneous subcellular Ca²⁺ oscillations (the criticality mechanism) (Nivala et al., 2012). SERCA2a pumping Ca²⁺ into SR and K⁺ channel repolarization of V_m return the Ca²⁺ and V_m domain cues toward their diastolic levels at which LCRs again begin to emerge, creating the ignition phase of the next AP cycle.

Self-Organized Criticality

Spatiotemporal self-organization across geometric scales (subcellular to cell-wide) is a manifestation of criticality that has been observed in excitable cells throughout nature (Stožer et al., 2019) including cultured astrocytes (Jung et al., 1998), immature oocytes (Lopez et al., 2012), and mouse cardiac ventricular myocytes (Nivala et al., 2012). Self-similar, scale-free distributions of parameters across wide scales that obey power law behavior (are ln–ln linear) are an indication of their self-ordered criticality (Bak, 1999).

It has recently been discovered that coupling of subcellular Ca²⁺ signals (cues) generated by the Ca²⁺ clock within isolated SAN cells to the cell surface membrane proteins during APs to elicit a change in V_m manifests long-range power law correlations (are self-similar) across species (Sirenko et al., 2021). Specifically, Ca²⁺ and V_m domain kinetic transitions (cues) during AP cycles in single SAN cells isolated from mice, guinea pigs, rabbits, and humans are self-similar to each other during APS and self-similar to trans-species APFIs of these cells in vitro, to heart rate in vivo, and to body mass (Sirenko et al., 2021).

Neurotransmitter stimulation of β-adrenergic receptor stimulation (βAR) or cholinergic receptor (CR)–initiated signaling modulates the AP firing rate and rhythm of SAN cells by impacting on coupled-clock protein functions, modulating the degree to which criticality (Ca²⁺ domain) and limit-cycle (V_m domain) mechanisms couple to each other during AP cycles (Maltsev and Lakatta, 2009; Lakatta et al., 2010). APFIs in rabbit SAN cells during autonomic stimulation vary over a fourfold range, from approximately 200 ms during βARs up to approximately 800 ms during CR stimulation (CRs) (Vinogradova et al., 2002; Lyashkov et al., 2009).

We hypothesized that transitions in V_m and Ca²⁺ domain cues during the diastolic AP ignition (Lyashkov et al., 2018) and recovery phases (Figure 1) of APs are (1) not only self-similar to each other in cells without autonomic receptor stimulation (control cells), but are self-similar to V_m and Ca²⁺ cues in other cells during CRs and during βARs and (2) that Ca²⁺ and V_m cues during APs are self-similar to APFI variabilities (and therefore self-similar to mean APFIs) regardless of the presence or absence or type of autonomic receptor stimulation. In other terms, we hypothesized that Ca²⁺ and V_m domain clock-coupling cues occurring during all APs are self-similar to each other, i.e., manifest long-range correlations in all isolated SAN cells within populations of cells that differ with respect to autonomic input, and that these Ca²⁺ and V_m cues during APs are also self-similar to the rate and rhythm of AP firing across the entire range of APFIs created by these cues in all isolated SAN cells.

FIGURE 1

Figure 1. Simultaneous recording of AP and Ca²⁺ in a 10 beat time series of APs (A) and the measured APFI and CaTFI (firing internal of AP-induced Ca²⁺ transient) (B); the definition of parameters measured in V_m and Ca²⁺ domains (C) and the self-similarity of the V_m and Ca²⁺ parameters during APs to each other and to the APFI (D).

To test these hypotheses, we studied a large population (n = 230) of single rabbit SAN cells to which we applied CRs [carbachol (CCh)] to one subset of cells, βARs [isoproterenol (ISO)] to another subset, and no autonomic receptor stimulation to a third subset of cells. This created three populations of SAN cells having APFIs distributed across the entire physiologic range. We measured intracellular Ca²⁺ or membrane potential in these cells to (1) characterize the times to ignition onset and times to 90% recovery of V_m and Ca²⁺ parameters during APs in AP time series and (2) to determine the correlations of these V_m and Ca²⁺ kinetic parameters to each other during APs, to APFI variability (APFIV; and therefore to mean APFIs). Thus, the data set to be analyzed consisted of 12 different kinetic parameters in each cell population (control, CCh, and ISO); six parameter means, three each in the Ca²⁺ and V_m domains; and six parameter variabilities (SDs) around the means. To determine the degree of self-similarity among V_m and Ca²⁺ domain parameters, we constructed density distribution plots and applied correlation, power law, and principal component (PC) analyses to Ca²⁺ and V_m domain data sets separately and to the combined Ca²⁺ and V_m data sets. We further extended our perspectives from cell population and single-cell levels downward to the molecular scale by performing numerical modeling simulation and analyzing variabilities of ion currents and Ca²⁺ with respect to APFI to determine whether these ion currents and Ca²⁺ also obeyed a power law across autonomic states.

Materials and Methods

The study was performed in accordance with the Guide for the Care and Use of Laboratory Animals published by the National Institutes of Health (NIH publication 85–23, revised 1996). The experimental protocols have been approved by the Animal Care and Use Committee of the NIH (protocol #457-LCS-2024). Materials and methods briefly presented here are detailed in Supplementary Material.

Isolation of Single Rabbit SAN Cells

Single, spindle-shaped, spontaneously beating SAN cells were isolated from the hearts of New Zealand rabbits (Charles River Laboratories, Wilmington, MA, United States) as described previously (Vinogradova et al., 2000).

Spontaneous APs Recordings

Time series of spontaneous APs were recorded in subsets of freshly isolated SAN cells using the perforated patch-clamp technique with Axopatch 200B patch-clamp amplifier (Axon Instruments) (Bogdanov et al., 2001) at 34°C ± 0.5°C. AP parameters (Figure 1) measured via a customized program (Lyashkov et al., 2018) were APFI, APD₉₀, and the time to ignition onset (TTIO) measured by the time at which diastolic membrane potential dV/dt accelerates to 0.15 V/s (Figure 1), which reflects the onset of the ignition phase of the AP cycle (Lyashkov et al., 2018).

Ca²⁺ Measurements

In another subset of SAN cells, AP-induced global CaTs and spontaneous LCRs (Figure 1) were measured at 34°C ± 0.5°C with a confocal microscope (Zeiss LSM510, Germany) in the line-scan mode (Vinogradova et al., 2004; Yang et al., 2012). The interval between the peaks of two adjacent AP-induced CaTs (Figure 1) is defined as CaT firing interval, which is highly correlated with the APFI, as demonstrated by simultaneous recordings of V_m and Ca²⁺ in a separate subset of cells (Figure 1). The LCR period is defined as the time from the peak of the prior AP-induced CaT to an LCR peak in diastole (Figure 1); the time to 90% decay of the CaT was defined as CaT₉₀.

Numerical Modeling

We performed numerical simulations using a modified Maltsev–Lakatta model that features the coupled-clock mechanism (Maltsev and Lakatta, 2009). The computer code for the original model is freely available and can be downloaded and run in CellML format¹ using the Cellular Open Resource software developed by Alan Garny at Oxford University in the United Kingdom (Garny et al., 2009) (for recent development of this software)². The original model could not be directly used for APFIV simulations because it is a system of first-order differential equations that is deterministic and showing no APFIV in limit-cycle oscillatory regime of AP steady firing. Thus, we modified the model to generate variability of AP waveforms by supplementing total membrane current (I_tot) with an additional randomly fluctuating current around its zero-mean value, known as perturbation current or I_per [as previously implemented by Henggui Zhang (Monfredi et al., 2014)]. Furthermore, we also performed an additional set of simulations with I_per added to Ca²⁺ release flux to mimic the effect of stochastic LCRs (Bogdanov et al., 2006; Monfredi et al., 2013). Using the resultant stochastic dynamical system of SAN cell, we simulated fluctuating APFI, ion currents, and Ca²⁺ dynamics for three conditions: (i) basal AP firing, (ii) during βARs with ISO (100 nM), and (iii) during CRs with CCh (100 nM). The effects of autonomic modulation (conditions ii and iii) were modeled as previously described (Maltsev and Lakatta, 2010), except modulation of I_CaL current by CCh that was modeled as described by Zaza et al. (1996). All model equations and parameters are provided in the Supplementary Material.

Experimental Design and Statistics

Supplementary Figure 1 illustrates schematic of the experimental design to assess long-range correlations of V_m and Ca²⁺ parameters during APs and APFI intervals in cells within and among populations of cells that differed with respect to autonomic input. V_m and Ca²⁺ parameter intervals (milliseconds) are presented as mean ± SD. APFIV within a time series is taken as standard deviation (SD) about the mean or as the coefficient of variation (CV) (the ratio of SD to the mean).

Analyses of V_m and Ca²⁺ parameter interval distributions measured in AP time-series in different cells determined the association among pairs of variables using Pearson’s correlations for both average and individual data (Howell, 2002). In cells in which Ca²⁺ was measured, the mean interval between AP-induced CaTs was usually longer than the mean APFI in cells in which APs were recorded (due to slight buffering effects of the fluorescent Ca²⁺ probe). To allow all the variables to be combined into a single analysis, we matched on the APFI variable Z scores in control, ISO, or CCh populations by “matchit” function in R (Ho et al., 2011).

Density estimates of within-cell standard deviations and means of parameters of cells within each autonomic state population are presented as non-parametric kernel estimates of probability density functions, scaled so that the total area under each curve is unity (Silverman, 1986).

In order to determine whether distributions of parameter means and SDs are self-similar, i.e., obeyed a power law suggesting fractal-like behavior, we constructed ln–ln plots of distributions of the AP and Ca²⁺ function means and SDs measured across the broad range of apparent steady states in the absence of, or the presence of, βARs or CRs (Kucera et al., 2000; Yaniv et al., 2013).

The relationships among all distributions of all parameters (means and SDs) were also assessed in principal component analyses (Johnson and Wichern, 2008).

In a few cells in which experimentally measured parameter interval distributions measured in the V_m and Ca²⁺ domains in the same cell and for numerical simulation of ion currents and Ca²⁺ prior to and during autonomic receptor stimulation Poincaré indices were employed to define long-range correlations among variables, a Poincaré plot graphs a parameter (n), in an AP time series on the x axis versus the same parameter of the succeeding AP (n + 1) on the y axis; i.e., one takes a sequence of parameters and plots each one against the following parameter (Huikuri et al., 2000).

When statistical inference was performed, p < 0.05 was considered statistically significant.

Results

Assessment of Self-Similarity of Ca²⁺ and V_m Kinetic Interval Parameters to Each Other During APs and to APFIV

Figure 1A illustrates a time series of APs in a SAN cell during which Ca²⁺ and V_m were simultaneously measured in the absence of autonomic receptor stimulation. Kinetic transitions in Ca²⁺ and V_m parameters as illustrated in Figure 1B were assessed during each AP. Figures 1C,D shows that V_m and Ca²⁺ parameters during AP time series are self-similar to each other and are also self-similar to the APFIs (r = 0.844 in this cell).

Figure 2 illustrates the time series of AP intervals in Figure 1A plotted as phase-plane diagrams in which V_m is depicted as a function of Ca²⁺ throughout each cycle. The times at which various channels are activated throughout the V_m/Ca loop are indicated. The times of ignition and 90% recovery are also indicated. The point resolution is 3.072 ms, and the point spread indicates the rates at which the electrochemical signal changes during each cycle. The dashed line in the figure marks the border between the disordered and ordered molecular activation. The arrows indicate the direction of the electrochemical signal emergence.

FIGURE 2

Figure 2. (A) The time series of AP intervals in Figure 1A plotted as phase-plane diagrams in which V_m is depicted as a function of Ca²⁺ throughout each cycle. The times at which various channels are activated throughout the V_m/Ca phase-plane loop are indicated. The times of ignition and 90% recovery are also indicated (B). The point resolution is 3.072 ms, and the spread indicates the rates at which the electrochemical signal changes during each cycle. The dashed line in the figure marks the border between the disordered and ordered molecular activation. The arrows indicate the direction of the electrochemical signal emergence.

A Poincaré plot (scatter graph) constructed from consecutive data points in a time series (Figure 3A) is a convenient tool that provides information on correlations (self-similarity) of data across the time series. The x axis defines the parameter (n) occurrence in milliseconds, and the y axis defines the parameter occurrence at (n + 1). The Poincaré plot in Figure 3A depicts the data of the time series of the cell in Figure 1. Note that although the means vary over a threefold range, all six parameter means (the Ca²⁺ and V_m interval parameters measured during APs and APFIs) in the absence of autonomic receptor stimulation are described by a line of identity, indicating their self-similarity across the AP time series. Quantitative analysis of short- and long-term variability in a given time series of observations entails fitting an ellipse to each cloud of data points within the Poincaré plot (Figure 3B): The length of a line describing the slope of the long axis of each ellipse is referred to as SD2 of the data points (c.f. Figure inset); the length of the line describing the slope of the short axis, which is perpendicular in direction to the long axis line, is referred to as SD1. Note in Figure 3C that the SD1 is self-similar to SD2 across the fourfold range of Ca²⁺ and V_m parameters. The center point of each ellipse, i.e., the intersection of SD1 and SD2, is the average interval between events (AP intervals or other parameters measured in the time series) within the time series. SD1/SD2 (Figure inset) informs on non-linear trends (unequal lengths of SD1 and SD2) across intervals within each ellipse.

FIGURE 3

Figure 3. Poincaré plots (A) and fitting of the Poincaré plot ellipse clouds (B), the relation of SD2s to SD1s (C) of the six parameters from simultaneously recorded during the 10 beats time series from Figure 1.

Figure 4 illustrates combined Poincaré plots of TTIO, APD₉₀, during APs and APFIs in time series of APs of two representative cells: one cell in control and during CRs by CCh and the other cell in control and during βARs stimulation by ISO.

FIGURE 4

Figure 4. Self-control AP recordings in two cells during CCh or ISO (A), Poincaré plots of the three measured V_m domain parameters in three autonomic states: CCh, control (two cells), and ISO (B, number of beats in each time series was 197, 397, 397, and 900, respectively, in carbachol, carbachol control, isoproterenol control, and isoproterenol). The inset shows an example of ellipse fitting.

Although the range of absolute values of kinetic interval parameters of cells depicted in the Poincaré plot in Figure 4B vary by 20-fold, all points (n = 5,673) are self-similar, i.e., fit by a single line (r = 0.992) with a slope of unity, passing nearly through the origin.

Supplementary Table 1 lists the SD1s, SD2s, SD1/SD2, the means of the TTIO and APD₉₀ intervals, and APFIs depicted in Figure 4B. Note also that the SD1s, SD2s, and SD1/SD2 of TTIO, APD₉₀, and APFIs progressively increase from ISO to control and markedly increase from control to CCh, creating degrees of non-linearity across the combined control, ISO, and CCh states, which is also reflected in the mean APFIs across the three states (Supplementary Table 1). The V_m transitions during APs across different autonomic states are self-similar to APFIs across these states (Figure 5A). Figure 5B shows the self-similarity of parameter means of SD1s of V_m parameters to their SD2s across autonomic states, indicating self-similarity of short-term (e.g., beat to beat) and long-term (e.g., rhythm across more than two beats) variabilities across autonomic states within the time series. Figure 5C shows the self-similarity of all parameter means depicted in Figures 4, 5A,B to their SDs across autonomic states.

FIGURE 5

Figure 5. Correlation between mean and SD of the three V_m parameters (A), correlation of Poincaré ellipse fitting SD1 and SD2 (B), and the self-similarity of APD90 and TTIO to APFI (C) across the three autonomic states, the same as in Figure 4.

The data in Figures 4, 5 demonstrate that over the entire range of physiologic APFIs from 192.7 ms in ISO to 305.5 ms in control, and to 910.3 ms in CCh, variabilities of TTIO and APD₉₀ measured during APs are self-similar to each other and are also self-similar to the variability of APFIs within the time series, and therefore self-similar to the mean APFI of the time series.

Self-similarity among V_m variables in the cells in Figure 4 across autonomic states in control and during ISO and CCh in Figures 2, 3A can also easily be ascertained from the shapes of their population density distributions (Figure 6). Note that in the cell superfused with CCh, the distributions of TTIO, APD₉₀, and APFIs (Figures 6A–C) are broader than in control or during ISO. Note also that the distributions of kinetic transitions during APs and APFI become more synchronized from CCh to ISO (Figure 6A–C). In other terms, the degree to which V_m and Ca²⁺ parameters are synchronized during APs increases from CCh to control to ISO, similar to the AP firing variabilities and mean APFIs.

FIGURE 6

Figure 6. Density distributions of the self-control V_m parameters, TTIO (A), APD90 (B) and APFI (C), measured in the same cell as in Figures 4, 5 (two control cells are combined). The density distributions are presented as non-parametric kernel estimates of probability density functions (Silverman, 1986), scaled so that the total area within each curve is unity.

Self-Similarity of Ca²⁺ and V_m Parameter Means to Each Other During APs and to APFI Variabilities and Mean APFIs Across Autonomic States in Different Cells

We next determined whether the self-similarity (long-range correlations) of Ca²⁺ to V_m parameters measured within the same cells as depicted in Figures 1–6 extends to populations of different cells within and among different autonomic states. To accomplish this, we applied CRs (CCh) to one subset of cells, βAR stimulation (ISO) to another subset, and no autonomic receptor stimulation to a third subset of cells. This created populations of SAN cells having APFIs distributed across the entire physiologic range.

Table 1 lists descriptive statistics (means and SDs) of the variabilities and means of Ca²⁺ and V_m kinetic parameter intervals measured in different cells within and among populations of cells that differ with respect to autonomic input. The mean SD of V_m and Ca²⁺ domain parameters listed in Table 1A gives the average time-series variability of each parameter among cells within each of the three cell populations (control, ISO, or CCh). The standard deviation of the SDs (SDSD) in Table 1A tells us how variable the SDs of each parameter are among cells within each cell population. The means of each parameter measured within a time series (Table 1B) tell us the average level of the parameter among cells within each of the three populations, and the SD of the means tells us the variability of the mean parameter levels among cells within each cell population.

TABLE 1

Table 1. (A) Mean of SDs and SD of SDs of AP and Ca²⁺ domain intervals among individual cells in each of the three steady state populations that differ with respect to autonomic receptor stimulation; (B) Mean ± SD of means of AP and Ca²⁺ domain interval in each cell population (A).

Figure 7A illustrates the distributions of SDs of parameters within Ca²⁺ and V_m domains during APs and of APFIs in cells listed in Table 1A of each of the three populations of cells that differed in autonomic state: control cells (n = 78), cells during superfusion with ISO (n = 27), and in cells superfused with CCh (n = 10). The self-similarity (long-range correlations) of the mean SDs of Ca²⁺ to V_m parameter transitions during ignition and recovery phases of APs across the wide range of APFIs induced by the type of autonomic receptor stimulation, or lack thereof, is evident in the self-similarity of their mean SD distribution shapes (Figure 7A).

FIGURE 7

Figure 7. Density distributions of selected parameter means (B) and SDs (A) of M and Ca²⁺ clock functions measured in different cells prior to and during autonomic receptor stimulation across the three groups of mean APFI steady states in Table 1. The density distributions are presented as non-parametric kernel estimates of probability density functions (Silverman, 1986), scaled so that the total area within each curve is unity. Both the mean and variability about the mean are concordant with each other across the three experimental groups in control and shift concordantly in response to autonomic receptor stimulation in all cells measured.

The distribution of the means listed in Table 1B is illustrated in Figure 7B. Note that the shapes of the distributions are self-similar to each other across the three different autonomic states. Note also that the shapes of the distribution of the means of a given parameter in Figure 7B are similar to the distribution of that parameter’s SDs in Figure 7A (because the interval distribution means stem from the distributions of their SDs). In other terms, the variability in the times at which parameters occur within a time series (their parameter SDs) determines what the mean interval of events in the time series will be.

Also note in Figure 7 that compared to cells not superfused with an autonomic receptor agonist (control cells) and those superfused with ISO, the shapes of the distributions of cells superfused with CCh are broad, indicating marked variability among CCh cells within the parameter distributions of both interval means and interval SDs.

Figure 8A illustrates the self-similarity of V_m parameter means and SDs to Ca²⁺ means and SDs across autonomic states. Heatmaps of the long-range correlations among V_m and Ca²⁺ parameter means are shown in Supplementary Figure 3. The long-range correlations (self-similarity) between the means of the means and means of SDs given in Table 1 is shown in Figure 8B.

FIGURE 8

Figure 8. (A) Plot of mean of means and SDs of V_m versus Ca²⁺ from Table 1. (B) Mean of means in Table 1 versus mean of the SDs in Table 1. (C) Means and SDs of V_m parameters versus means and SDs of corresponding Ca²⁺ parameters for all 230 cells. (D) Means and SDs of ignition and recovery versus means and SDs of cycle length for all 230 cells (In panel (A), M.M = mean of means; M.SD = mean of SDs; CL = cycle length; Ign = Ignition; Rec = recovery).

Self-Similarity of Ca²⁺ and V_m Parameters Among All Individual Cells Within and Among the Three Autonomic States

Two-by-two correlations of V_m and Ca²⁺ parameter means and SDs of all 230 cells that comprised the three different cell populations in Figure 7 and Table 1 are highly significant (Table 2).

TABLE 2

Table 2. Correlation matrices of V_m or Ca²⁺ parameters measured in different cells within three different autonomic states (summary data listed in Table 1).

Figure 9 shows ln–ln plots of the distributions of means and SDs of Ca²⁺ and V_m parameters of all cells across the three autonomic states. The piecewise linear fit of the data in each panel is largely driven by the data from cells superfused with CCh that manifested broad interval distributions and high mean APFIs in Figures 4–7. Figure 8C shows that the means and SDs of V_m parameters measured during APs in individual cells (n = 115) are self-similar to Ca²⁺ means and SDs measured during AP in other individual cells (n = 115) across the three autonomic states.

FIGURE 9

Figure 9. Distributions of SDs and means over all groups illustrating self-similarity (fractal-like behavior). The regression lines are fit as a single piecewise linear model with a join point at the center of the interval with the highest frequency.

Although, as noted above, the V_m and Ca²⁺ parameters within and among cells during CCh superfusion were more broadly distributed than those during control or during ISO, the correlations between V_m and Ca²⁺ parameters among all cells within the CCh superfused population of cells were extremely strong for most parameters (Supplementary Figure 3). Weaker but still significant correlations between times to 90% recovery and other variables are observed in CCh superfused cells (likely because times to 90% are the most difficult parameters in the data set to measure accurately).

Figure 8D shows that V_m and Ca²⁺ parameters measured during APs are self-similar to APFI means and SDs across the three autonomic states.

Correlation of Ca²⁺ and V_m Domain Parameter Means in Individual Cells to Their SDs Within and Across Autonomic States

The relationship between mean AP firing rate and its SD is known to be non-linear (Monfredi et al., 2014). Whereas the relationships of all Ca²⁺ and V_m parameter means relative to their corresponding SDs measured in the combined set of data derived from different populations of cells in control or during ISO or CCh superfusion are non-linear (Figures 10A–C), the ln–ln plots of these combined data (Figure 10D–F) are linear, indicating their self-similarity across all 230 cells that differed by autonomic state.

FIGURE 10

Figure 10. Means versus SDs for all cells on a linear scale for (A) V_m, (B) Ca²⁺, and (C) V_m and Ca²⁺, and using a logarithmic scale for (D) V_m, (E) Ca, and (F) V_m and Ca²⁺.

Principal Component Analyses

Next, we employed principal component analyses to determine whether the self-similarity of parametric measures within the entire data set of variables could be summarized by a smaller set of principal components that contain most of the information in all the variables. PCs are linear combinations of the original variables, and each PC is statistically independent of the others: the first PC explains as much of the total variability in the data as possible, the second PC as much of the remaining variability, and so on. Highly self-similar parameters within the complete data set are explained by the sum of the first few PCs.

In a PC analysis of SDs of the six variables (three in the V_m domain and three in the Ca²⁺ domain), the first three PCs explained 91.4% of the total variation within the entire SD data set (Supplementary Table 2 and Figure 11A). Similarly, in a PC analysis of the means of all six measured variables means (three in the V_m domain and three in the Ca²⁺ domain) in a PC analysis, the first two PCs explained 92.8% of the variability in the data set of all six means (Supplementary Table 2 and Figure 11B). Finally, in a PC analysis of all 12 variables (six means and six SDs), the first three PCs explained 88% of the variability within the total (means plus SDs) data set (Supplementary Table 2 and Figure 11C).

FIGURE 11

Figure 11. Plots of principal components: (A) SDs only, (B) means only, (C) means and SDs.

Because a smaller set of PCs can explain a substantial proportion of the total variability in each set of V_m or Ca²⁺ domain means, SDs, and means and SDs, means that these distributions of V_m and Ca²⁺ parameters measured during an AP and APFIs in control cells and different cells superfused with ISO or CCh are each self-similar to each other. In other terms, Ca²⁺ and V_m domain functions operative within the SAN cells coupled-clock system manifest self-similar scale-free characteristics, i.e., kinetic fractals of each other, across the entire physiologic range of APFIs.

Numerical Model Simulations of APFIV, Major Ion Currents, and Ca²⁺

Variability of V_m and Ca²⁺ parameters measured experimentally in cells within and across autonomic states is linked to the respective variabilities of clock molecular availability to respond to V_m and Ca²⁺ cues (Figure 2) that cannot be directly measured experimentally during AP firing. To gain further insight into the variability of these biophysical mechanisms, we performed numerical modeling simulations. The APFIV was generated by SAN cell model described as a stochastic dynamical system, i.e., a dynamical system [deterministic Maltsev–Lakatta model (Maltsev and Lakatta, 2009)] subjected to the effects of noise current, I_per (see section “Materials and Methods” and Supplementary Material for details). I_per amplitude was tuned for the model APFIV to match that measured experimentally under respective experimental conditions. We investigated two scenarios of noise generation: when I_per was added to I_tot or when I_per was added to Ca²⁺ release flux current in each of the three autonomic states: (i) basal AP firing, (ii) ISO 100 nM, and (iii) CCh (100 nM). Variability of six major currents was simulated and analyzed: I_f, I_NCX, I_Kr, I_CaL, I_CaT, and I_KACh. Variability of [Ca] under cell membrane was also simulated during the three autonomic states. For all items, we measured variability of their peak amplitudes and amplitudes at −40 mV during DD.

Model simulation results are presented in Figure 12, with their numerical values given in Supplementary Tables 6, 7. Regardless of the type of noise generation (via Ca²⁺ or I_tot), it affected the variability of ion currents and Ca²⁺ the same way, and the predicted variabilities for many parameters differed substantially from that of APFI:

(1) I_f variability was substantial: in the basal state and in ISO I_f variability was similar to or larger than APFIV; the variability of I_f decreased in CCh.

(2) I_NCX variability was also substantial: at −40 mV, it was substantially larger than that of APFIV (except in CCh when I_per was added to I_tot); variability of I_NCX peak amplitude (negative, during AP upstroke) was similar to that of APFI in the basal state and ISO, but became reduced in CCh.

(3) I_Kr variability was substantially less than APFIV under all conditions.

(4) I_CaT variability was the largest among ion currents, being similar to that of I_NCX, in CCh when I_per was added to Ca²⁺ release.

(5) Peak I_CaL variability was always less than that of APFIV; at −40 mV, it was greater in ISO than in basal state and CCh.

(6) Variability of Ca²⁺ release flux in the basal state and ISO at −40 mV was greater than or similar to APFIV; in CCh, variability of Ca²⁺ release flux was less than that of APFI.

FIGURE 12

Figure 12. Results of analysis of coefficient of variation (CV) of major ion currents I_f, I_Kr, I_CaL, I_CaT, I_KACh, and Ca²⁺ simulated by Maltsev–Lakatta coupled-clock model with noise current (I_per) added to either total current I_tot (A–C) or Ca release flux (D–F).

Some components exhibited power law behavior over a wide range of APFI over all conditions tested (Figure 13).

FIGURE 13

Figure 13. I_NCX (A), I_CaT (B), I_f (C), and [Ca]_sub (D) exhibited power law behavior (a linear dependence in the ln–ln plot here) over a wide range of APFI over all conditions tested. Noise (I_per) was added to I_tot in Maltsev–Lakatta coupled-clock model. Similar dependencies were found when I_per was added to Ca release flux (not shown).

We next determined whether self-similarity across autonomic states observed for experimental data during the ignition phase is also applied to simulated ion currents or Ca²⁺ data during this time of the cycle, i.e., at −40 mV. To this end, we applied the statistical tests utilized for experimental data to simulated data (for the scenario when I_per was added to I_tot).

Simulated ion currents and Ca²⁺ amplitudes during AP ignition (−40 mV) across the three autonomic states are self-similar to each other, strongly correlated to each other, as were experimentally measured parameters (Table 2). These two-by-two correlations of all the simulated components are listed in Supplementary Table 3. Selected examples of these correlations are shown in Figure 14. As I_NCX is fully determined by V_m and Ca²⁺, at a fixed voltage (−40 mV) it is fully determined by only Ca²⁺. That is why we have 100% correlation of I_NCX and Ca²⁺. I_CaT strongly correlated with Ca²⁺ variations, because the stronger Ca²⁺ signal is linked to the higher DD rate and hence stronger (time-dependent) activation of I_CaT. Surprisingly, I_KACh amplitude was also highly correlated with Ca variations.

FIGURE 14

Figure 14. The correlation between I_NCX (A), I_CaT (B), I_f (C), I_Kr (D), or I_KACh (E) versus [Ca]_Sub at –40 mV. Note that I_KACh presents only during CCh.

Variations in I_f and I_Kr at −40 mV did not depend on variations of Ca²⁺, but their mean values strongly depended on Ca²⁺ across the autonomic states. I_f activation and I_Kr deactivation are early DD mechanisms and do not seem to interplay with Ca²⁺ at the ignition onset at −40 mV in a given cycle.

Figure 15A illustrates the Poincaré plots of many simulated parameters (APFI and TTIO; I_NCX, I_CaT, I_Kr, [Ca], I_CaL, I_KACh, and I_f, all at −40 mV). Although the ranges of absolute values of these simulated parameters substantially vary, all simulated parameters are self-similar, i.e., fit by a single line (r = 0.998) with a slope of unity, resembling Poincaré relationship of experimentally measured AP parameters (Figures 1, 4).

FIGURE 15

Figure 15. Poincaré plots of the simulated I_NCX, I_CaT, I_CaL, I_Kr, I_KACh, I_f, and [Ca]_Sub at –40 mV, APFI and TTIO from one model cell during ISO, CCh or at basal condition in ln–ln plot (A), and the correlation between the means and SDs of these parameters (B). Note that I_KACh only presents during carbachol.

Figure 15B illustrates ln–ln plots of the relationships of the means of simulated components to their SDs. Note that this relationship follows power law behavior just as did the relationship of experimentally measured means for AP parameters versus their SDs (Figure 10).

Finally, in PC analyses of the simulated parameters in Figure 15B, the first two PCs accounted for 94% of the variation in the eight variables (Supplementary Figure 4).

Discussion

We measured membrane potential and Ca²⁺ times to onsets of AP ignition during diastole and times to 90% recovery during APs and APFIs in three populations of single, isolated rabbit SAN cells that differed with respect to autonomic input: those in which CRs were stimulated by CCh, those in which βARs were stimulated with ISO, and in untreated (control) cells. Absolute values of times to AP ignition onsets and to 90% recovery intervals in Ca²⁺ and V_m domains during APs and APFIs differed within and among individual cells in cell populations and different autonomic states and differed markedly among cells of the three populations of cells with differential autonomic receptor stimulation. Our novel finding is that although differing markedly in absolute values, Ca²⁺ and V_m parameters were self-similar to each other during APs and self-similar to APFIs, not only within and among different cells within each of the three populations of cells studied, but remarkably, among all cells, regardless of the autonomic receptor stimulation profile. Thus, Ca²⁺ and V_m domain kinetic transitions (intervals) during APs, individual APFI, and mean APFIs within AP time series manifest long-range correlations (self-similar scale-free correlations, i.e., obey power law) across the entire broad range of APFIs, regardless of whether autonomic receptors of these cells are stimulated and regardless of the type of autonomic receptor stimulation.

The degree to which molecular activation states within each clock and between clocks are synchronized during APs determines when the next AP will occur, i.e., the APFIV and mean APFIs within a cell and across the entire population of single SAN cells studied: the higher the degree of order (self-organized activation of clock molecules), the more ordered and less variable the aggregate of kinetic functions, the least variability of APFIs, and the shorter the mean APFI; vice versa, the lower the degree of order among clock molecular activation states, the lower the aggregate synchronization among clock molecular functions, the greater the variability of APFIs, and the longer the mean AP cycle interval.

Self-similar or fractal-like beating rate variability among cardiac cells in culture has been previously identified in a number of studies but only when cells were confluent or electrically connected to each other. This behavior has been attributed to influences of tonic or phasic resetting of membrane potential or to mechanical factors via cell-to-cell connections (Clay and DeHaan, 1979; Jongsma et al., 1983; Kucera et al., 2000). Our novel observation is self-similarity of V_m and Ca²⁺ domain intervals during APs and APFIs across diverse populations of single SAN cells that were not physically connected to each other.

Thus, self-similar distributions of order that have been demonstrated to occur in other instances throughout nature (Bak, 1999) also exist within SAN cell coupled-clock system functions. We interpret this power law behavior of SAN cell functions to result from concordant gradations of self-organized order (synchronization) of clock molecular activation across the entire physiologic range of APFIs.

Clock Molecular Activation Cues

Voltage, time, Ca²⁺, cAMP signaling, and PKA and CaMKII-dependent clock protein phosphorylation are the cues that regulate the activation kinetics of molecules that control pacemaker functions in single SAN cells (Supplementary Figure 2; Lakatta et al., 2003, 2006, 2008, 2010; Maltsev and Lakatta, 2008; Yaniv et al., 2015). Some coupled-clock system proteins are activated by Ca²⁺, e.g., SERCA2; others by V_m and cAMP binding, e.g., HCN channels (DiFrancesco and Tortora, 1991) and other cyclic nucleotide-regulated channels; or by Ca²⁺ and V_m, e.g., NCX, or by phosphorylation and Ca²⁺, e.g., phospholamban and ryanodine receptors and AC type 8, whereas the activation states of still other coupled-clock system proteins are modulated by V_m, Ca²⁺, and phosphorylation, e.g., L-type and some K⁺ channels.

Both voltage and Ca²⁺ activation cues oscillate in amplitude throughout each AP cycle and command rapid responses from clock molecules. The degree to which activation status of molecules of a given species is synchronized at any given time following the prior AP determines the ensemble response of that molecular species to its activation cues. It is well documented that following a synchronizing event, e.g., the occurrence of an AP, activation states of molecules underlie AP cycle transition through variably inactivated states, altering the availability to respond to a subsequent activation cue. Our new concept of synchronization of functional cues is based on the idea that the coupled-clock system inheres (inevitably) some degree of disorder that stems from its key constituent proteins operating (stochastically switching) intrinsically within their conformational flexibilities and heterogeneity. The balance of order/disorder is linked to molecule interactions (i.e., effectiveness of their respective cues) that allow them to operate cooperatively as an ensemble or system with various degrees of synchronization (i.e., order) that is reflected in respective variability of the output function of the system, i.e., APFIV in our case.

Thus, we interpret the experimentally measured concordant behavior of surface membrane and Ca²⁺ regulatory functions during AP cycles across the entire physiologic range of AP cycles to reflect a concordance in the degrees of activation of molecules that drive these regulatory functions. Importantly, Ca²⁺ and V_m cues not only regulate the synchronization clock of molecular activation states but are also regulated by the degree of synchronized activation of molecules determined by these cues (recursion). Because membrane and Ca²⁺ clocks become coupled in the context of the electrochemical signal that waxes and wanes to cause the AP cycle, the extent of self-organized molecular activation within each clock indirectly affects self-organization of molecular activation of the other clock operating within the coupled-clock system. And because scaling of mean APFIs among all cells is self-similar to APFIV among cells, AP firing variability and mean APFI are determinants of the Ca²⁺ and V_m cues that determine kinetic intervals during an AP: a recursive, feed-forward process.

Mean APFI and APFIV Are Not Only Regulated by but Also Regulate the Degree to Which Clock Molecular Functions Are Synchronized

Changes in Ca²⁺ and V_m cues during an AP determine not only the characteristics of that AP but also when the next AP will occur and the mean APFI within an AP time series.

A prolongation of the mean APFIs, itself, contributes to the concurrent increase in the APFIV at a long mean APFI: because an increase in mean APFI reduces net Ca²⁺ influx and indirectly reduces Ca^2+/CaMKII-AC–dependent phosphorylation of Ca²⁺ cycling proteins, reducing the SR Ca²⁺ cycling kinetics and increasing the variability of LCR periods.

Characteristics of the AP that are determined by availability of M clock molecules to respond to a change in membrane potential both directly and indirectly entrain the Ca²⁺ and M clock activation: as the mean AP interval shortens, less time elapses between APs, and therefore at shorter intervals, less time is required than at longer intervals for molecules to retain (remember) the synchronizing influences imparted by the preceding AP. This causes the relationship of mean APFI to APFIV of isolated SAN cells to be non-linear (Figure 5 and Supplementary Figure 2), as originally demonstrated by Zaza and Lombardi (2001) and later by Monfredi et al. (2014). Conversely, as time following a prior AP increases, the effectiveness of the Ca²⁺ activation cue, itself, wanes because the cell Ca²⁺ level and SR Ca²⁺ load become reduced because of time-dependent Ca²⁺ efflux from the cell. We may speculate therefore that during long AP cycles, fewer molecules of some molecular species are available to respond to Ca²⁺ activation cues.

Gradations of self-organized molecular activation within and between clocks regulate the APFI rhythm, i.e., the APFIV. In other terms, the average APFI, kinetics of the AP, AP-triggered Ca²⁺-transient, LCR periods and DD kinetics, and beat-to-beat variability of these parameters measured in the present study are readouts of the relative extents to which of clock molecules become activated and the degree to which the clocks are coupled. When the degree to which Ca²⁺ and M clocks kinetics are coupled or synchronized is low, the AP firing rate is slow, and APFIV is high, e.g., during CRs. Conversely, when the degree of coupling or synchronization of the Ca²⁺ and V_m kinetics of the two clocks is high, e.g., during βAR stimulation, AP firing is rapid, and APFIV is low.

So, What Factors Affect the Degree of Synchronization of Clock Molecules?

Concordant degrees of self-similar synchronization of M and Ca²⁺ clock kinetic functions reflect concordant gradations of activation states of specific molecules that govern these functions and how these cues change throughout an AP cycle.

AP Firing Rate and Rhythm Synchronization of Clock Molecules

The AP that emerges from the diastolic ignition events is, itself, the most potent integrator or synchronizer, not only of surface membrane electrogenic molecules, but also of Ca²⁺ clock functions: a synchronized global cytosolic CaT that ensues following synchronous activation of voltage-dependent L-type Ca²⁺ channels is created by synchronized Ca²⁺-induced, Ca²⁺ release from SR via ryanodine receptor activation (Song et al., 1998; Wang et al., 2001; Lakatta, 2004; Zhou et al., 2009).

The efficacy of V_m and Ca²⁺ activation cues that oscillate as electrochemical signal that underlies the V_m change during AP cycle varies with the AP cycle interval or period: shorter periods (i.e., faster AP firing rates or shorter APFIs) are more effective than longer periods (i.e., slower AP firing rates or longer APFIs), because during very long AP cycles, Ca²⁺ activation states of some molecules become more unsynchronized. At very short times following a large voltage oscillation (i.e., an AP), many molecules of a given molecular species in relatively inactivated state may not optimally respond to activation cues (e.g., impaired excitability/non-excitability). As the time following a prior activation increases, although a subpopulation of molecules of given species may regain full ability to respond to activation cues, substantial variability in the activation status of other molecules of that species still may exist, limiting the number of molecules that can respond to (be recruited by) an activation cue. Our results provide novel clues to the cellular basis for the observation that an AP occurrence, itself, influences the range of APFIs that immediately follow it (Nolasco and Dahlen, 1968). The AP, itself, indirectly affects all Ca²⁺ clock functions because it regulates net cell Ca²⁺ balance. Functions of M clock molecules that underlie the generation of an AP indirectly regulate the availability for SR Ca²⁺ cycling by modulation of the level of cell Ca²⁺, the SR “oscillatory substrate.” Thus, M clock functions also indirectly regulate LCR periods and sizes via their impact on the “steady state” intracellular Ca²⁺ level. When the average interval between APs becomes prolonged, a reduction net Ca²⁺ influx into efflux from the cell (Lakatta, 2004) reduces the cytosolic [Ca²⁺], the rate of Ca²⁺ pumping into SR, and the SR Ca²⁺ load. These reductions, in turn, prolong the average time from the prior AP occurrence for spontaneous local diastolic ryanodine receptor activation to occur within SAN cell local microdomains; the randomness of spontaneous local diastolic ryanodine receptor activation occurring within these microdomains also increases, broadening the distribution of LCR periods and shifting these to longer times at a long AP cycle (Figure 4).

Thus, the degree of variability in activation states of M and Ca²⁺ clock molecules that emerges over time following their synchronization by the prior AP is implicated in the cycle length dependence of variability of Ca²⁺ and M clock functions measured here (Table 1 and Figures 3–5). Heartbeat variability in vivo and APFIV of isolated SAN cells in vitro indicate that neither autonomic input to SAN cells, nor functions intrinsic to the SAN cell coupled-clock system, respectively, achieve a steady state from one beat to the next.

Ca²⁺-Dependent Synchronization of Clock Molecules

The local [Ca²⁺], itself, also serves as a powerful synchronizer of clock molecular function: ordered/disordered Ca²⁺ regulation has been recently reported for ryanodine receptor–mediated Ca²⁺ releases (Maltsev et al., 2019).

Studies in permeabilized SAN cells, in which Ca²⁺ clock function is preserved, but M clock function is abolished, and therefore APs cannot occur and do not influence LCR periodicity, clearly demonstrate that in a fixed, physiologic, free [Ca²⁺], LCR occurrences are random when the free [Ca²⁺] is low and that LCR periodicity emerges as the free [Ca²⁺] in the system is increased due in part to an increase in the Ca²⁺ charge of the SR capacitor (Sirenko et al., 2013). The intracellular concentration of the oscillatory substrate, Ca²⁺, itself is regulated, in part, by the SAN cell transmembrane Na⁺ gradient and membrane potential (Sirenko et al., 2013, 2016).

cAMP Activation or Phosphorylation of Clock Proteins Modulates the Synchronization of and Response to Activation Cues

Autonomic receptor stimulation modulates both the activation cues and responses of clock molecules to these cues. The impact of autonomic receptor signaling on the effectiveness of clock coupling occurs over several AP cycles and is reflected in time-dependent transitions in the AP firing rate and rhythm. The kinetics and stoichiometry of increases in PKA activity in response to gradations in βAR stimulation predict the kinetics and stoichiometry of concurrent time-dependent increases in AP firing rate (Yaniv et al., 2014b). Prior studies (Lyashkov et al., 2009; Yang et al., 2012) have demonstrated that gradations in the phosphorylation status of phospholamban at Ser¹⁶ across the three autonomic state mean APFIs of cell populations in the present study (Supplementary Figure 3) strikingly resemble gradations of the means of APFIs and APFIVs observed across these autonomic states measured in the present study (Table 1).

The extent to which clock molecules respond to Ca²⁺ or V_m activation cues during an AP is modulated by βAR- or CR-dependent phosphorylation of many of the same proteins that drive SAN cell automaticity in the absence of autonomic receptor activation. These βARs or CRs impact on the memory of the extent to which clock molecule activation had been synchronized during the prior AP. βAR or CR modulation has two facets: (1) a direct effect, due to cAMP or phosphorylation-dependent activation of clock proteins and (2) an indirect effect by altering the APFI, which alters the cell Ca²⁺ activation cues that are directly modulated by autonomic receptor stimulation. Specifically, βARs and CRs, respectively, not only reduce or increase mean APFI, but also, respectively, shift variability within distributions of Ca²⁺ and V_m functions in the same direction (Sirenko et al., 2016) as the shift in mean APFI.

Numerical Modeling

Because we experimentally measured characteristics of APs in populations of single cells that differed by autonomic input status, we were able to glean insights not only into APFIV in an “average” cell, but also into populations of cells isolated from SAN tissue. Embracing SAN function at the cell population level resonates with recent studies of SAN function at the tissue level (Bychkov et al., 2020; Fenske et al., 2020) that have revealed a novel understanding of the SAN impulse as an emergent property created by a collective of variable interactions among heterogeneous cell populations within the SAN tissue. On the other end, APFIV per se also emerges at smaller, subcellular scales, due to variability in the functions of individual molecules, such as ion channels, transporters, and pumps, individually and in complex cooperation with each other. These molecular functions cannot be measured directly in our single-cell experiments. Thus, we employed numerical modeling to extend our perspectives from cell populations and single-cell levels downward to the molecular scale. Such broader consideration of variability makes sense when we approach pacemaker function as a multiscale phenomenon (Qu et al., 2011) featuring free scale and fractal-like characteristics (Weiss and Qu, 2020). Considering the SAN cell as a stochastic dynamic system, we examined variability of major ion currents and submembrane [Ca²⁺] during different autonomic states that created a broad range of APFIs, which were measured experimentally.

Our simulations indicate that the APFIV of some ion currents and submembrane Ca²⁺ can be close to that of the APFI itself, but also can be substantially lower or higher than the APFIV, depending on the presence or absence of autonomic receptor stimulation, and the time during the AP cycle (Figure 12): components such as I_f, I_NCX, I_CaT, and Ca²⁺ exhibit substantial cycle to cycle variability, whereas I_CaL, I_KACh, and I_Kr show less or moderate variability. This behavior reflects complex non-linear recursive interactions of V_m and Ca²⁺ that couple the clocks that drive the system (Lyashkov et al., 2018) and as such cannot be directly and definitively interpreted in cause-effect terms. Nevertheless, our simulations confirmed ion channel behavior that could have been envisioned. For example, independent of the nature of the noise added (to I_tot or to Ca²⁺ release flux): components contributing to DD (I_f, I_NCX, and I_CaT) exhibit a larger cycle to cycle variability, whereas components contributing to generation of all-or-none AP characteristics exhibit less variability (peak I_CaL and I_Kr). This result is in accord with the idea of order/disorder transitions during AP cycle (Figure 2, order/disorder dash line), i.e., the DD manifests disorder and transition to order, and hence, a larger variability, and then following the AP upstroke, the AP itself manifests order and hence less variability.

However, our simulations also provided some unexpected, interesting results: I_NCX and submembrane Ca²⁺ amplitudes during DD followed a power law relationship over a wide range of APFI under all conditions tested, indicating that self-similar scale-free or fractal-like behavior is operative within the coupled-clock mechanism via I_NCX (Figure 13). There is also likely to be a secondary effect of I_NCX amplitude itself on DD acceleration in a recursive fashion. Indeed, an increase in I_NCX is expected to accelerate DD but, at the same time, is further accelerated by the very same acceleration it imparts to the ignition mechanism (i.e., diastolic I_CaL and Ca²⁺-induced Ca²⁺ release). This self-acceleration ignition mechanism results in the power law behavior predicted by the model.

Peak I_f also followed a power function over a broad range of APFIs, albeit in a noisy manner manifesting some extremely long APFIs in CCh. The current that fluctuated most with respect to variabilities in APFI turned out to be I_CaT, the peak amplitude of which also reflects DD dynamics. When DD is rapidly accelerating (as it does at shorter cycles), I_CaT peak quickly activates and achieves a higher peak amplitude, and vice versa, when DD is slower at longer APFIs, I_CaT becomes inactivated over time without achieving a peak. In log–log plots, this forms a straight line for almost the entire range of APFIs (power law behavior), except extremely long APFIs when variations in (already) slow DD dynamics have almost no effect on I_CaT peak. Thus, simulations of biophysical components within the pacemaker cell system exhibited a power law behavior over a wide range of APFIs that encompasses the broad range of APFIs measured experimentally.

Observation of model simulations through analytic lenses applied to experimental data indicated that (as was the case for the experiment data) the simulated variables are self-similar to each other across a broad range of APFIs within the three autonomic states (Figure 15, Supplementary Figure 4, and Supplementary Table 3).

Model Limitation

Our I_CaL model was adopted from the Kurata model (Kurata et al., 2002), and although it does include Ca²⁺-dependent inactivation mechanism, it lacks I_CaL facilitation described in 2000 by Mangoni et al. (2000). Also, our numerical model of I_f lacks dynamic regulation by cAMP (DiFrancesco, 1999), Ca²⁺ activated K⁺ channels and store-operated channels that may potentially contribute to APFIV.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Ethics Statement

The animal study was reviewed and approved by the Animal Care and Use Committee of the National Institutes of Health (protocol #457-LCS-2024).

Author Contributions

EL and VM designed the project. DY, AL, IZ, YY, TV, and BZ performed the experiments and analyzed the results. DY, CM, and ST did the statistical analyses. DY and VM performed the numerical simulation. DY and EL wrote the manuscript with the support from YY and VM. All authors contributed to the article and approved the submitted version.

Funding

This work was supported partially by the Intramural Research Program of the NIH, National Institute on Aging, and by the ISF [No. 330/19].

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors wish to acknowledge the assistance of Dr. Jia-Hua Qu for assistance with the PC analysis illustrations, and sincerely appreciate Loretta Lakatta, Robert Monticone, Tracy Oppel, and Ruth Sadler for their editorial assistance. This manuscript has been released as a pre-print at bioRxiv (https://www.biorxiv.org/content/10.1101/2020.09.01.277756v2) (Yang et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.612770/full#supplementary-material

Abbreviations

AP, action potential; APD₉₀, action potential duration from AP overshoot to 90% repolarization; APFI, AP firing interval; APFIV, variability of AP firing intervals; TTIO, time from the previous AP overshoot to the ignition onset when dV/dt = 0.15 (V/S); CaTFI, firing internal of AP-induced Ca²⁺ transient; CaT₉₀, 90% decay time of AP-induced Ca²⁺ transient; LCR, local Ca²⁺ releases; SR, sarcoplasmic reticulum; SD, standard deviation; CV, coefficient of variation; PC, principal component; I_CaL, high voltage-activated, L-type Ca²⁺ current; I_f, hyperpolarization-activated funny current; I_KACh, acetylcholine-activated K⁺ current; I_Kr, K⁺ current exhibiting strong inward rectification; I_NCX, Na⁺-Ca²⁺ exchanger (NCX) current; β ARs, β-adrenergic receptor stimulation; CRs, cholinergic receptor stimulation; CCh, carbachol; ISO, isoproterenol; HR, heart rate; SAN, sinoatrial node.

Footnotes

1. ^ http://models.cellml.org/workspace/maltsev_2009
2. ^ http://www.opencor.ws/

References

Bak, P. (1999). How Nature Works: The Science of Self-Organized Criticality. New York, NY: Springer-Verlag.

Google Scholar

Bogdanov, K. Y., Maltsev, V. A., Vinogradova, T. M., Lyashkov, A. E., Spurgeon, H. A., Stern, M. D., et al. (2006). Membrane potential fluctuations resulting from submembrane Ca2+ releases in rabbit sinoatrial nodal cells impart an exponential phase to the late diastolic depolarization that controls their chronotropic state. Circ. Res. 99, 979–987. doi: 10.1161/01.RES.0000247933.66532.0b

CrossRef Full Text | Google Scholar

Bogdanov, K. Y., Vinogradova, T. M., and Lakatta, E. G. (2001). Sinoatrial nodal cell ryanodine receptor and Na+-Ca2+ exchanger: molecular partners in pacemaker regulation. Circ. Res. 88, 1254–1258. doi: 10.1161/hh1201.092095

PubMed Abstract | CrossRef Full Text | Google Scholar

Bychkov, R., Juhaszova, M., Tsutsui, K., Coletta, C., Stern, M. D., Maltsev, V. A., et al. (2020). Synchronized cardiac impulses emerge from multi-scale, heterogeneous local calcium signals within and among cells of heart pacemaker tissue. JACC Clin. Electrophysiol. 6, 907–931. doi: 10.1016/j.jacep.2020.06.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Clay, J. R., and DeHaan, R. L. (1979). Fluctuations in interbeat interval in rhythmic heart-cell clusters. Role of membrane voltage noise. Biophys. J. 28, 377–389. doi: 10.1016/S0006-3495(79)85187-5

CrossRef Full Text | Google Scholar

DiFrancesco, D. (1999). Dual allosteric modulation of pacemaker (f) channels by cAMP and voltage in rabbit SA node. J. Physiol. Pharmacol. 515, 367–376. doi: 10.1111/j.1469-7793.1999.367ac.x

PubMed Abstract | CrossRef Full Text | Google Scholar

DiFrancesco, D., and Tortora, P. (1991). Direct activation of cardiac pacemaker channels by intracellular cyclic AMP. Nature 351, 145–147. doi: 10.1038/351145a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenske, S., Hennis, K., Rotzer, R., Brox, V. F., Becirovic, E., Scharr, A., et al. (2020). cAMP-dependent regulation of HCN4 controls the tonic entrainment process in sinoatrial node pacemaker cells. Nat. Commun. 11:5555. doi: 10.1038/s41467-020-19304-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Garny, A., Noble, D., Hunter, P. J., and Kohl, P. (2009). CELLULAR OPEN RESOURCE (COR): current status and future directions. Philos. Trans. A Math. Phys. Eng. Sci. 367, 1885–1905. doi: 10.1098/rsta.2008.0289

PubMed Abstract | CrossRef Full Text | Google Scholar

Ho, D., Imai, K., King, G., and Stuart, E. (2011). MatchIt: nonparametric preprocessing for parametric causal inference. J. Stat. Softw. 42, 1–28.

Google Scholar

Howell, D. C. (2002). Statistical Methods for Psychology, 5th Edn. Duxbury, MA: Thomson Learning.

Google Scholar

Huikuri, H. V., Makikallio, T. H., Peng, C. K., Goldberger, A. L., Hintze, U., and Moller, M. (2000). Fractal correlation properties of R-R interval dynamics and mortality in patients with depressed left ventricular function after an acute myocardial infarction. Circulation 101, 47–53. doi: 10.1161/01.cir.101.1.47

CrossRef Full Text | Google Scholar

Johnson, R. A., and Wichern, D. W. (2008). Applied Multivariate Statistical Analysis, 6th Edn. Upper Saddle River, NJ: Pearson.

Google Scholar

Jongsma, H. J., Tsjernina, L., and de Bruijne, J. (1983). The establishment of regular beating in populations of pacemaker heart cells. A study with tissue-cultured rat heart cells. J. Mol. Cell. Cardiol. 15, 123–133. doi: 10.1016/0022-2828(83)90288-2

CrossRef Full Text | Google Scholar

Jung, P., Cornell-Bell, A., Madden, K. S., and Moss, F. (1998). Noise-induced spiral waves in astrocyte syncytia show evidence of self-organized criticality. J. Neurophysiol. 79, 1098–1102. doi: 10.1152/jn.1998.79.2.1098

PubMed Abstract | CrossRef Full Text | Google Scholar

Kucera, J. P., Heuschkel, M. O., Renaud, P., and Rohr, S. (2000). Power-law behavior of beat-rate variability in monolayer cultures of neonatal rat ventricular myocytes. Circ. Res. 86, 1140–1145. doi: 10.1161/01.res.86.11.1140

CrossRef Full Text | Google Scholar

Kurata, Y., Hisatome, I., Imanishi, S., and Shibamoto, T. (2002). Dynamical description of sinoatrial node pacemaking: improved mathematical model for primary pacemaker cell. Am. J. Physiol. 283, H2074–H2101.

Google Scholar

Lakatta, E. G. (2004). Beyond Bowditch: the convergence of cardiac chronotropy and inotropy. Cell Calcium 35, 629–642. doi: 10.1016/j.ceca.2004.01.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakatta, E. G. (2021). Heartbeat music. Heart Rhythm 18, 811–812. doi: 10.1016/j.hrthm.2021.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakatta, E. G., Maltsev, V. A., Bogdanov, K. Y., Stern, M. D., and Vinogradova, T. M. (2003). Cyclic variation of intracellular calcium: a critical factor for cardiac pacemaker cell dominance. Circ. Res. 92, e45–e50.

Google Scholar

Lakatta, E. G., Maltsev, V. A., and Vinogradova, 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. doi: 10.1161/CIRCRESAHA.109.206078

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakatta, E. G., Vinogradova, T., Lyashkov, A., Sirenko, S., Zhu, W., Ruknudin, A., et al. (2006). The integration of spontaneous intracellular Ca2+ cycling and surface membrane ion channel activation entrains normal automaticity in cells of the heart’s pacemaker. Ann. N. Y. Acad. Sci. 1080, 178–206. doi: 10.1196/annals.1380.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Lakatta, E. G., Vinogradova, T. M., and Maltsev, V. A. (2008). The missing link in the mystery of normal automaticity of cardiac pacemaker cells. Ann. N. Y. Acad. Sci. 1123, 41–57. doi: 10.1196/annals.1420.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopez, L., Piegari, E., Sigaut, L., and Dawson, S. P. (2012). Intracellular calcium signals display an avalanche-like behavior over multiple lengthscales. Front. Physiol. 3:350. doi: 10.3389/fphys.2012.00350

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyashkov, A. E., Behar, J., Lakatta, E. G., Yaniv, Y., and Maltsev, V. A. (2018). Positive feedback mechanisms among local Ca Releases, NCX, and ICaL ignite pacemaker action potentials. Biophys. J. 114, 1176–1189. doi: 10.1016/j.bpj.2017.12.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyashkov, A. E., Vinogradova, T. M., Zahanich, I., Li, Y., Younes, A., Nuss, H. B., et al. (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. doi: 10.1152/ajpheart.01340.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, A. V., Maltsev, V. A., Mikheev, M., Maltseva, L. A., Sirenko, S. G., Lakatta, E. G., et al. (2011). Synchronization of stochastic Ca2+ release units creates a rhythmic Ca2+ clock in cardiac pacemaker cells. Biophys. J. 100, 271–283. doi: 10.1016/j.bpj.2010.11.081

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, A. V., Stern, M. D., and Maltsev, V. A. (2019). Mechanisms of calcium leak from cardiac sarcoplasmic reticulum revealed by statistical mechanics. Biophys. J. 116, 2212–2223. doi: 10.1016/j.bpj.2019.04.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, V. A., and Lakatta, E. G. (2008). Dynamic interactions of an intracellular Ca2+ clock and membrane ion channel clock underlie robust initiation and regulation of cardiac pacemaker function. Cardiovasc. Res. 77, 274–284. doi: 10.1093/cvr/cvm058

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, V. A., and Lakatta, E. G. (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. doi: 10.1152/ajpheart.01118.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltsev, V. A., and Lakatta, E. G. (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. doi: 10.1152/ajpheart.00783.2009

PubMed Abstract | CrossRef Full Text | Google Scholar

Mangoni, M. E., Fontanaud, P., Noble, P. J., Noble, D., Benkemoun, H., Nargeot, J., et al. (2000). Facilitation of the L-type calcium current in rabbit sino-atrial cells: effect on cardiac automaticity. Cardiovasc. Res. 48, 375–376. doi: 10.1016/s0008-6363(00)00182-6

CrossRef Full Text | Google Scholar

Monfredi, O., Lyashkov, A. E., Johnsen, A. B., Inada, S., Schneider, H., Wang, R., et al. (2014). Biophysical characterization of the underappreciated and important relationship between heart rate variability and heart rate. Hypertension 64, 1334–1343. doi: 10.1161/HYPERTENSIONAHA.114.03782

PubMed Abstract | CrossRef Full Text | Google Scholar

Monfredi, O., Maltseva, L. A., Spurgeon, H. A., Boyett, M. R., Lakatta, E. G., and Maltsev, V. A. (2013). Beat-to-beat variation in periodicity of local calcium releases contributes to intrinsic variations of spontaneous cycle length in isolated single sinoatrial node cells. PLoS One 8:e67247. doi: 10.1371/journal.pone.0067247

PubMed Abstract | CrossRef Full Text | Google Scholar

Nivala, M., Ko, C. Y., Nivala, M., Weiss, J. N., and Qu, Z. (2012). Criticality in intracellular calcium signaling in cardiac myocytes. Biophys. J. 102, 2433–2442. doi: 10.1016/j.bpj.2012.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Nolasco, J. B., and Dahlen, R. W. (1968). A graphic method for the study of alternation in cardiac action potentials. J. Appl. Physiol. 25, 191–196. doi: 10.1152/jappl.1968.25.2.191

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, Z., Garfinkel, A., Weiss, J. N., and Nivala, M. (2011). Multi-scale modeling in biology: how to bridge the gaps between scales? Prog. Biophys. Mol. Biol. 107, 21–31. doi: 10.1016/j.pbiomolbio.2011.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Shivkumar, K., Ajijola, O. A., Anand, I., Armour, J. A., Chen, P. S., Esler, M., et al. (2016). Clinical neurocardiology defining the value of neuroscience-based cardiovascular therapeutics. J. Physiol. 594, 3911–3954. doi: 10.1113/JP271870

PubMed Abstract | CrossRef Full Text | Google Scholar

Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. London: Chapman & Hall/CRC.

Google Scholar

Sirenko, S., Maltsev, V. A., Yaniv, Y., Bychkov, R., Yaeger, D., Vinogradova, T. M., et al. (2016). Electrochemical Na+ and Ca2+ Gradients drive coupled-clock regulation of automaticity of isolated rabbit sinoatrial nodal pacemaker cells. Am. J. Physiol. Heart Circ. Physiol. 311, H251–H267. doi: 10.1152/ajpheart.00667.2015

PubMed Abstract | CrossRef Full Text | Google Scholar

Sirenko, S., Yang, D., Li, Y., Lyashkov, A. E., Lukyanenko, Y. O., Lakatta, E. G., et al. (2013). Ca2+-dependent phosphorylation of Ca2+ cycling proteins generates robust rhythmic local Ca2+ releases in cardiac pacemaker cells. Sci. Signal. 6:ra6. doi: 10.1126/scisignal.2003391

PubMed Abstract | CrossRef Full Text | Google Scholar

Sirenko, S. T., Tsutsui, T., Tarasov, K. V., Yang, D., Wirth, A. N., Maltsev, V. A., et al. (2021). Self-similar synchronization of calcium and membrane potential transitions during action potential cycles predict heart rate across species. JACC Clin. Electrophysiol. doi: 10.1016/j.jacep.2021.02.016 [Epub ahead of print].

CrossRef Full Text | PubMed Abstract | Google Scholar

Song, L. S., Sham, J. S., Stern, M. D., Lakatta, E. G., and Cheng, H. (1998). Direct measurement of SR release flux by tracking ‘Ca2+ spikes’ in rat cardiac myocytes. J. Physiol. 512(Pt 3), 677–691. doi: 10.1111/j.1469-7793.1998.677bd.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Stožer, A., Markovic, R., Dolensek, J., Perc, M., Marhl, M., Rupnik, M. S., et al. (2019). Heterogeneity and delayed activation as hallmarks of self-organization and criticality in excitable tissue. Front. Physiol. 10:869. doi: 10.3389/fphys.2019.00869

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinogradova, T. M., Bogdanov, K. Y., and Lakatta, E. G. (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. doi: 10.1161/hh0102.102271

PubMed Abstract | CrossRef Full Text | Google Scholar

Vinogradova, T. M., Bogdanov, K. Y., Yang, D., Kuschel, M., Cheng, H., and Xiao, R. P. (2000). Sinoatrial node pacemaker activity requires Ca2+/calmodulin-dependent protein Kinase II activation. Circ. Res. 87, 760–767. doi: 10.1161/01.res.87.9.760

CrossRef Full Text | Google Scholar

Vinogradova, T. M., Zhou, Y. Y., Maltsev, V., Lyashkov, A., Stern, M., and Lakatta, E. G. (2004). Rhythmic ryanodine receptor Ca2+ releases during diastolic depolarization of sinoatrial pacemaker cells do not require membrane depolarization. Circ. Res. 94, 802–809. doi: 10.1161/01.RES.0000122045.55331.0F

CrossRef Full Text | Google Scholar

Wang, S. Q., Song, L. S., Lakatta, E. G., and Cheng, H. (2001). Ca2+ signalling between single L-type Ca2+ channels and ryanodine receptors in heart cells. Nature 410, 592–596. doi: 10.1038/35069083

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, J. N., and Qu, Z. (2020). The sinus node: still mysterious after all these years. JACC Clin. Electrophysiol. 6, 1841–1843.

Google Scholar

Yang, D., Lyashkov, A. E., Li, Y., Ziman, B. D., and Lakatta, E. G. (2012). RGS2 overexpression or Gi inhibition rescues the impaired PKA signaling and slow AP firing of cultured adult rabbit pacemaker cells. J. Mol. Cell. Cardiol. 53, 687–694. doi: 10.1016/j.yjmcc.2012.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, D., Morrell, C. H., Lyashkov, A. E., Tagirova, S., Zahanich, I., Yaniv, Y., et al. (2020). Ca²⁺ and membrane potential transitions during action potentials are self-similar to each other and to variability of ap firing intervals across the broad physiologic range of AP intervals during autonomic receptor stimulation. bioRxiv [Preprint]. doi: 10.1101/2020.09.01.277756

CrossRef Full Text | Google Scholar

Yaniv, Y., Ahmet, I., Liu, J., Lyashkov, A. E., Guiriba, T. R., Okamoto, Y., et al. (2014a). Synchronization of sinoatrial node pacemaker cell clocks and its autonomic modulation impart complexity to heart beating intervals. Heart Rhythm 11, 1210–1219. doi: 10.1016/j.hrthm.2014.03.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaniv, Y., Ganesan, A., Yang, D., Ziman, B. D., Lyashkov, A. E., Levchenko, A., et al. (2015). Real-time relationship between PKA biochemical signal network dynamics and increased action potential firing rate in heart pacemaker cells: kinetics of PKA activation in heart pacemaker cells. J. Mol. Cell. Cardiol. 86, 168–178. doi: 10.1016/j.yjmcc.2015.07.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Yaniv, Y., Lyashkov, A. E., and Lakatta, E. G. (2013). The fractal-like complexity of heart rate variability beyond neurotransmitters and autonomic receptors: signaling intrinsic to sinoatrial node pacemaker cells. Cardiovasc. Pharm. Open Access. 2:111.

Google Scholar

Yaniv, Y., Lyashkov, A. E., Sirenko, S., Okamoto, Y., Guiriba, T. R., Ziman, B. D., et al. (2014b). Stochasticity intrinsic to coupled-clock mechanisms underlies beat-to-beat variability of spontaneous action potential firing in sinoatrial node pacemaker cells. J. Mol. Cell Cardiol. 77, 1–10. doi: 10.1016/j.yjmcc.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaza, A., and Lombardi, F. (2001). Autonomic indexes based on the analysis of heart rate variability: a view from the sinus node. Cardiovasc. Res. 50, 434–442. doi: 10.1016/s0008-6363(01)00240-1

CrossRef Full Text | Google Scholar

Zaza, A., Robinson, R. B., and Difrancesco, D. (1996). Basal responses of the L-type Ca2+ and hyperpolarization-activated currents to autonomic agonists in the rabbit sino-atrial node. J. Physiol. Pharmacol. 491, 347–355. doi: 10.1113/jphysiol.1996.sp021220

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., Zhao, Y. T., Guo, Y. B., Xu, S. M., Bai, S. H., Lakatta, E. G., et al. (2009). Beta-adrenergic signaling accelerates and synchronizes cardiac ryanodine receptor response to a single L-type Ca2+ channel. Proc. Natl. Acad. Sci. U.S.A. 106, 18028–18033. doi: 10.1073/pnas.0906560106

PubMed Abstract | CrossRef Full Text | Google Scholar