Self-similar action potential cycle-to-cycle variability of Ca2+ and current oscillators in cardiac pacemaker cells

Variability of heart pacemaker cell action potential (AP) firing intervals (APFI) means that pacemaker mechanisms do not achieve equilibrium during AP firing. We tested whether mechanisms that underlie APFI, in rabbit sinoatrial cells are self-similar within and across the physiologic range of APFIs effected by autonomic receptor stimulation. Principal Component Analyses demonstrated that means and variabilities of APFIs and local Ca2+ releases kinetics, of AP induced Ca2+-transient decay times, of diastolic membrane depolarization rates, of AP repolarization times, of simulated ion current amplitudes, are self-similar across the broad range of APFIs (264 to 786 ms). Further, distributions of both mean APFIs and mean Ca2+ and membrane potential dependent coupled-clock function kinetics manifested similar power law behaviors across the physiologic range of mean APFIs. Thus, self-similar variability of clock functions intrinsic to heart pacemaker cells determines both the mean APFI and its interval variability, and vice versa.


Introduction
The heart's beating rate and rhythm are regulated by the rate and rhythm of spontaneous action potential (AP) firing of sinoatrial nodal pacemaker cells (SANC) residing within sinoatrial node (SAN) tissue. Although, by convention, the number of heart beats per minute is casually referred to as "the" in vivo heart rate (HR), analyses of successive R-R intervals in an EKG timeseries detects substantial beat-to-beat interval variability, referred as to HR variability, largely attributable to autonomic input to the SANC.
Evidence subsequently emerged to indicate that spontaneous AP firing in single, isolated SANC is controlled by a mutually entrained, coupled-oscillator system: (1) an ensemble of surface membrane electrogenic molecules or "M clock" that directly controls the membrane potential and trans-membrane ion flux, and indirectly regulates intracellular Ca 2+ cycling; and (2) an intracellular "Ca 2+ clock", the sarcoplasmic reticulum (SR) and its decorator proteins, that directly control intracellular Ca 2+ cycling, and indirectly regulate transmembrane ion flux (Lakatta, Maltsev, & Vinogradova, 2010;. Subsequent studies suggested that APFIV of single SANC originates, in part, within the Ca 2+ clock (Carre et al., 1993;Monfredi et al., 2013;Zaniboni et al., 2014). Specifically, it has been demonstrated that APFIV within a given "steady state" emanates from the extent to which stochastic local Ca 2+ releases (LCR) generated by the Ca 2+ clock become spatiotemporally synchronized during diastolic intervals (Monfredi et al., 2013), and to the fidelity or effectiveness of clock coupling characteristics in a given apparent steady state (Lyashkov et al., 2018;Tsutsui et al., 2018;Yaniv, Lyashkov, et al., 2014).
We hypothesized the idea that APFIV of single, isolated SANC emerges from concordant cycle-to-cycle variability in the extent to which molecules operating within the coupled oscillator system become activated (due to synchronization of their activation states) during an AP cycle.
We further hypothesized that mean APFIs are self-similar across the entire physiologic range. We applied autonomic receptor stimulation, i.e., β adrenergic receptor stimulation (bARs) or cholinergic receptor stimulation (CRs), in order to effect a broad physiological range of mean APFIs. Although the concordance in cycle-to-cycle variability of activation of clock molecules cannot be directly measured when SANC are synchronously firing APs, concordance of the functional parameters of surface membrane voltage and intracellular Ca 2+ , that indirectly inform on the concordance of activation states of their molecular underpinnings, can be measured during spontaneous AP firing,. Such functional mechanisms include periodic spontaneous diastolic LCRs, kinetics of AP induced Ca 2+ -transient decay, diastolic membrane depolarization and AP repolarization.
We recorded the mean and cycle-to-cycle variability of APFI and the mean and variability of the aforementioned M and Ca 2+ clock functions during AP firing in single freshly isolated and cultured adult rabbit SANC (Yang, Lyashkov, Li, Ziman, & Lakatta, 2012). We also utilized a numerical coupled-clock model (V. A. Maltsev & Lakatta, 2009) to simulate cycle-to-cycle variability of selected ion current amplitudes by voltage clamping the model using experimentally measured AP time-series.

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 number. 85-23, revised 1996). The experimental protocols have been approved by the Animal Care and Use Committee of the National Institutes of Health (protocol #034LCS2016).

Single, isolated rabbit SANC Isolation and culture.
Single, spindle-shaped, spontaneously beating SANC were isolated from the hearts of New Zealand rabbits (Charles River Laboratories, Wilmington, MA) as described previously (Vinogradova et al., 2008). We also employed cultured SANC (c-SANC), in which phosphorylation becomes reduced compared to freshly isolated SANC (f-SANC), and mean APFI in c-SANC increases to about twice that of f-SANC and remains stable for several days (Yang et al., 2012).

Spontaneous AP Recordings.
Spontaneous APs (number of beats >512) were recorded in subset-sets of SANC using the perforated patch-clamp technique with an Axopatch 200B patchclamp amplifier (Axon Instruments) (Bogdanov et al., 2006) at 340.5 o C. AP parameters measured via a customized program (Bogdanov et al., 2006) were APFI, APD90, and the time from maximum diastolic polarization to the onset non-linear DD (TNLDD), which reflects the onset of the ignition phase of the AP cycle (Lyashkov et al., 2018).

Ca 2+ Measurements
. Subsets of f-and c-SANC were loaded with the Ca 2+ indicator fluo-4/AM (Thermo Scientific) (Vinogradova et al., 2004;Yang et al., 2012). AP initiated global Ca 2+ transients and spontaneous LCR were measured at 340.5 o C with a confocal microscope (Zeiss LSM510, Germany) in the line-scan mode. Ca 2+ -Image processing, data analysis, and presentation were performed using our original programs written in IDL software  (Interactive Data Language, Harris Corporation). The interval between the peaks of two adjacent AP-triggered Ca 2+ -transients is defined as Ca 2+ -transient firing interval (CaTFI). This approach was validated by the extremely high correlation between APFI and CaTFI measured simultaneously in a separate subset of cells in which both Ca 2+ and membrane potential were simultaneously measured (Fig. 1).The LCR period is defined as the time from the peak of the prior AP-induced Ca 2+ transient to an LCR peak in diastole.

Numerical modeling.
We developed a new "in silico AP clamp" method to predict individual realistic ion current time-series featuring beat-to-beat variability. In essence, we dynamically voltage-clamped numerical SANC models with respective experimental AP recordings. It is important to distinguish this new "in silico AP clamp" method from the traditional experimental AP clamp method that applied clamps in an individual real cell to previously measured AP waveforms in the presence of a specific ion channel blocker published previously (Zaza, Micheletti, Brioschi, & Rocchetti, 1997). In traditional AP clamp methods, the electronic current passing the cell membrane in this condition compensates the inhibited current to generate the same AP shape and has been interpreted as the specific current during real AP firing. While the elegant traditional method permits assessment of physiologically relevant ion current profiles, it is rarely used because of its technical difficulty and limitations, linked to poor selectivity of ion channel blockers (and their use-, and voltage-dependent action) and also voltage divergence (or escape) of real membrane voltage from AP clamp command, especially during fast potential changes (i.e. the action potential upstroke). The escape happens mainly because of series resistance intrinsically present in whole cell (or perforated) patch clamp mode. Our new "in silico AP clamp" method lacks these limitations, because we clamp the virtual membrane instantly (with a time step of 0.01ms) to experimentally measured AP waveform in the set of differential equations of the model, and the model generates all individual ion currents without any pharmacological blockers.
A "basal state firing model" of rabbit SANC (V. A. Maltsev & Lakatta, 2009) (http://models.cellml.org/workspace/maltsev_2009), which was subsequently modified to simulate the effect of CRs (V. A. ) based on the graded CRs levels (Lyashkov et al., 2009) was employed to simulate ion currents and their beat-to-beat variability during AP firing. We performed these simulations under two conditions: prior to and in the presence of CRs with carbachol (CCh, 100 nM). The maximum SR Ca 2+ pumping rate and maximum conductance for funny current (If), L-type Ca 2+ current (ICaL), and acetylcholine-activated K + current (IKACh) under these conditions were tuned in the model (V. A. .

Statistics.
A firing interval variable is produced in each data set, APFI and CaTFI. APFI and AP waveform parameters, CaTFI and LCR periods are presented as Mean ± SE. Beat-to-beat variability is taken as standard deviation (SD) about the mean, or as the coefficient of variation (CV, the ratio of SD to the mean). Parameter density estimates are nonparametric kernel estimates of probability density functions, scaled so that the total area under each curve is unity (Silverman, 1986). A Poincaré Plot graphs a parameter (n) on the x-axis versus the same parameter (n+1, the parameter of the succeeding AP) on the y-axis, i.e. one takes a sequence of parameters and plots each one against the following parameter (Huikuri et al., 2000) .
Statistical analyses of experimentally measured parameters were performed both within and among data subsets in which the AP or the Ca 2+ parameters were separately measured. Specific analyses were conducted to determine: the association among pairs of variables using Spearman's correlations for average data, and Pearson's correlation for both average and individual data (Howell, 2002); the relationships among all the variables (means and SDs), using principal components analyses (Johnson & Wichern, 2008), for the AP or the Ca 2+ data sets separately, and also for the combined AP and Ca 2+ parameter data sets. Because the mean CaTFI in experiments in which Ca 2+ was measured was usually longer than APFI in cells in which APs were recorded, to allow all the variables to be combined in a single principal components analysis, we matched on the APFI variable Z-scores within each 5 groups by matchit function in R (Ho, Imai, King, & Stuart, 2011). We also constructed Ln-Ln plots of distributions of the AP and Ca function means measured across the broad range of apparent steady states in the absence and presence of bARs or CRs to determine whether this distribution are self-similar, i.e., confirmed to Power-Law, suggesting fractal-like behavior (Kucera, Heuschkel, Renaud, & Rohr, 2000;Y. Yaniv, A. E. Lyashkov, & E. G. Lakatta, 2013a), A p value<0.05 was considered statistically significant, when paired t-test or ANOVA was applied.

Simultaneous AP and Ca 2+ signals recorded in the same cell. Prior to experiments in
which AP and Ca 2+ were measured separately in different cells comprising a large cell population, we performed simultaneous recordings of AP and Ca 2+ in a small subset of cells at baseline and following autonomic stimulation. Fig. 1A illustrates an example of AP's and Ca 2+ -transients measured simultaneously in a single SANC. Fig. 1B shows that APFI and CaTFI in cells in which both Ca 2+ and membrane potential were simultaneously measured are highly corelated.
AP measurements prior to and in response to autonomic receptor stimulation in the same cell. We also measured AP parameters prior to and in response to autonomic receptor stimulation in the same cell. Fig. 2 illustrates the concordance of means of AP parameters in the same cell prior to and in response to bARs or CRs. Note the concordant shifts in means of the AP parameters from control to autonomic receptor stimulation. Poincaré Plots (Fig. 2G~I) illustrate concordant beat-to-beat variability within and among different measured parameters in both c-and f-SANC measured in the same cell prior to and during bARs or CRs. Note the similar trends between in both fresh and cultured cells.

Distributions of Ca 2+ or M clock functions measured in control and during autonomic receptor stimulation in different cells. Table 1 lists means and SDs of AP parameters in cells in
which membrane potential was recorded, and means and SDs of parameters in other cells loaded with the Ca 2+ probe, in each of the 5 apparent "steady states" in control and in different cells during autonomic receptor stimulation. Note that the data set encompasses a broad range of APFI measured in 348 cells from 264ms (in f-SANC during bARs) to 786ms (in f-SANC during the steady response to CRs). and concordantly shift to smaller times (or higher frequencies, i.e. kinetics of the measured functions became accelerated) vs. control; and the distribution density hovering about the shorter mean APFI about the means narrows, (i.e. a higher dominant mean frequency, resulting from an increased number of events occurring at that frequency). Conversely, in response to CRs, all parameters (means and SDs) shift to longer times than in control, exhibit the lowest mean frequency and highest cycle-to-cycle variability recorded within the entire data set. Statistical analyses of the average data in Table 1 indicate a strong association not only between the means and SDs of a given parameter across the 5 experimental groups, but also among all parameters across all groups. Across the groups, the means and SDs tend to increase together.

Statistical tests for coherency among distributions of M and Ca 2+ clock kinetic functions
across the entire data set. To statistically evaluate the remarkable degree of concordance among means and SDs of coupled clock kinetic parameters among different cells in control, and concurrent concordant shifts of these parameters as the mean APFI shifts in response to autonomic stimulation (Table 1, Figs. 3-5), we looked deeper into individual cells and performed simple correlations among the means of parameters listed in Table 1 (Table 2). Almost all functions (all means and SDs), albeit measured in different cells prior to and during bARs or CRs, were highly correlated with each other in all the cases. Note that this correlation holds not only for functions among different cells within the same (AP or Ca 2+ ) data subset, but also for functions measured among different cells in the other (Ca 2+ or AP) data subset (Table 2).
When means and SDs for all measured functions in both AP and Ca 2+ data subsets were included in the principal components analyses model, the first 3 principal components accounted for 96.1% or 89.6% of total variation among means and SDs respectively (Table 3), indicating strong concordance among all variables within the data set. When both means and SDs of the AP and Ca 2+ data sets were included in the model, the first 3 principal components accounted for 86.8% of the total variance among the 12 variables listed in Table 1 (Table 3).

Means and SDs of different functions measured in different subsets of cells across a broad physiologic range of mean APFI are self-similar (obey a Similar Power-law).
We determined whether the remarkable, graded coherence observed among means and distributions of clock functional parameters within and among groups in which mean APFI's differ in Tables 1-3 and Figs.3-5 might indicate that the relationships of both mean APFIs and mean clock functions and their SDs are self-similar within each steady state and across the 5 "steady states" in Table 1.
The Ln-Ln transform plots of distributions for all cells for means and SDs listed in Table 1, in control and in different cells during bARs or CRs indicate that all measured parameters within the data set are self-similar and exhibit bi-fractal-like Power law behavior. Interestingly, the range of APFI about the bi-fractal-like intersection (Fig. 6, Ln5.7 or 300msec and above) is similar to that at which mean APFI and mean APFIV relationship begins to transition from linear to nonlinear ( Fig. 6A) and in Fig. 5C, which depicts a more detailed APFI and APFIV relationship across a broad range of autonomic input shown in Fig. 5A. The slopes of the Ln-Ln distributions plots for 6 means within Ca 2+ and AP data subsets at APFIs above an approximate mean cycle length of about 300ms (Ln=5.7, frequency of 3.34Hz) are negative and tightly grouped (Fig. 6G). This frequency range broadly encompasses f-SANC during stimulation of muscarinic receptors by the cholinergic agonist CCh and c-SANC in control (  (Fig. 5C).

Numerical model simulations of cycle-to-cycle variability of selected ion currents.
In addition to measuring and analyzing mean and variability of AP functions, we explored the cycleto-cycle variability of the ion currents that underlie respective APs and APFIV. Although it is not possible to measure ion currents experimentally during spontaneous AP firing, their characteristics can be predicted by numerical modelling that stimulates currents under AP clamp using respective experimentally measured AP time-series (see Methods).
We examined experimental AP-time series recordings of ten f-SANC in the basal state (control) and in the presence of CRs with CCh (100 nM), i.e. when variability of ion currents is expected to be maximal on the basis of experimentally measured cycle variability (Fig.5 & Table   1). Fig. 7A shows the experimentally measured AP waveform in a representative cell in control and in response to CCh. The cycle-to-cycle simulated current peak amplitudes of selected currents, i.e., ICaL, the inwardly rectifying K + current (IKr), If, and Na + -Ca 2+ exchanger (NCX) current (INCX) in each AP cycle during spontaneous AP firing are shown in panel B. IKACh is not shown in control as it's not activated (Lyashkov et al., 2009). Diastolic INCX was measured at a membrane potential of -45mV (Fig. 7C) , i.e. near ICaL activation threshold . This specific voltage (and therefore timing) at which INCX is measured is crucial in assessing clock coupling, because at any time during the AP cycle, INCX depends on both membrane potential and submembrane [Ca 2+ ] values. INCX measured at a fixed (i.e. instant) membrane potential level of -45 mV indirectly, but explicitly, informs on the diastolic [Ca 2+ ] beneath the cell membrane and thus the ensemble diastolic Ca 2+ release signal at the AP ignition onset (Lyashkov et al., 2018). As stated above (Table 1), experimentally measured APFI variability was markedly enhanced by CCh (Panel D two left bars). Cycle-to-cycle variability in the amplitudes of both INCX and If also markedly increased to a similar extent as the experimentally measured cycle-to-cycle APFI variability measured in the same cells in which the current amplitudes were simulated. In contrast, simulated cycle-to-cycle variability of ICaL and IKr remained almost unchanged. Average numerical values of mean APFI measured experimentally and simulated currents in the same cells in response to CCh are listed in Table 4.

Discussion
Our results demonstrate (1) that all measured coupled clock functional parameters of single, isolated SANC measured during an AP cycle exhibit concordant cycle-to-cycle variabilities within a given apparent steady state (i.e. manifest concordant degrees of order), and (2) that cycle-tocycle variability of all parameters shifts concurrently and concordantly across a broad range of steady state mean APFI, evoked by perturbations of autonomic receptor signaling (Table 1, Figs. 3-5). Specifically, bARs and CRs not only respectively, reduce or increase mean APFI but also APFIV shifted in the same direction. We interpreted these results to indicate that bARs increases and CRs reduces the order within and among these coupled-clock functions that underlie the AP cycle. In other terms, bARs increases and CRs decreases, respectively, the degree of synchronization of pacemaker cell clock regulatory functions about shorter (bARs) or longer (CRs) mean APFI.

Concordant beat-to-beat variability of M and Ca 2+ clock regulatory functions during an AP cycle reflects concordant gradations of activation states of specific molecules that govern
these functions. Clock molecular protein activation during each AP cycle arises from numerous internal and external cues that include voltage, time, Ca 2+ , cAMP signaling, and PKA and CaMKIIdependent clock protein phosphorylation (Lakatta, Maltsev, Bogdanov, Stern, & Vinogradova, 2003;Lakatta et al., 2010;Lakatta et al., 2006;Lakatta, Vinogradova, & Maltsev, 2008;Yaniv et al., 2015). The kinetics of state transitions of molecules that regulate pacemaker functions in single SANC are governed by a series of rate constants that are modulated by Ca 2+ , phosphorylation and voltage (Fig. 8). Both voltage and Ca 2+ activation cues command rapid responses from clock molecules, and these cues oscillate in amplitude throughout each AP cycle. The degree to which state activation of molecules of a given species is synchronized at any given time following a 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 that underlie AP cycle transition through variably inactivated states, altering the availability of given molecules to respond to a subsequent activation cue. Our new concept of synchronization of the 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 allows them to operate cooperatively as an ensemble or system with various degree of synchronization (i.e. order) that is reflected in respective variability of the output function of the system, i.e. APFIV in our case.
The protein cooperation and synchronization are executed via dependencies of activation transitions on numerous common factors, including Ca and membrane potential. Thus, the AP that originates from the aforementioned diastolic ignition events is, itself, the most potent integrator or synchronizer of clock molecule functions: not only of surface membrane electrogenic molecules; but also of Ca 2+ -induced Ca 2+ release ryanodine receptor activation, resulting in a synchronized global cytosolic Ca 2+ transient. that ensues following synchronous activation of voltage-dependent L-type Ca 2+ channels (Lakatta, 2004;Song, Sham, Stern, Lakatta, & Cheng, 1998;Wang, Song, Lakatta, & Cheng, 2001;Zhou et al., 2009). Ca 2+ itself also serves as a powerful synchronizer: order/disorder regulation has been recently reported for ryanodine receptor -mediated Ca 2+ releases (A. V. Maltsev, Stern, & Maltsev, 2019).
The influence of the mean APFI on APFIV becomes reduced as the mean AP interval shortens, because at short intervals less time elapses between APs, and therefore less time is required to retain (remember) the synchronizing influences on the "factors" impacted by the preceding AP than at longer AP intervals. This causes the relationship of mean APFI to APFIV of isolated SANC to be non-linear (Fig. 5), as originally demonstrated by ZaZa (Zaza & Lombardi, 2001). Thus, gradations of self-organized molecular activation within and between clocks, not only regulate the mean APFI on a beat-to-beat basis, but also the APFIV, i.e. APFI rhythm. In other terms, the average APFI, kinetics of AP, AP-triggered Ca 2+ -transient, LCR periods and diastolic depolarization kinetics, and beat-to-beat variability of these parameters measured here are readouts of the relative extents to which of clock molecules become activated.
At very short times following a large voltage oscillation (during 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 sub-population of molecules of given species may regains full ability to respond to activation cues, substantial variability in activation states of other molecules of that species still may exist, limiting the number of molecules that can respond to (be recruited by) an activation cue.
The intracellular concentration of the oscillatory substrate, Ca 2+ , itself is regulated, in part, by the SANC transmembrane Na + gradient and membrane potential Sirenko et al., 2013). As time following a prior AP increases, the effectiveness of the Ca 2+ activation cue, itself, wanes, because the cell Ca 2+ level and SR Ca 2+ load become reduced, due to time dependent Ca 2+ 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 2+ activation cues.
Thus, the efficacy of membrane potential and Ca activations cues that oscillate during electrochemical gradients that underlie an AP cycle varies with the AP cycle period: shorter periods (i.e. faster AP firing rates or shorter AP cycle lengths) are more effective than longer periods (i.e. slower AP firing rates or longer AP cycle lengths), because during very long AP cycles, Ca 2+ activation states of some molecules become more unsynchronized.

Important role of clock protein phosphorylation on self-organization of both molecular activation cues and activation states of molecules that regulates APFIV and mean APFI. Studies
in permeabilized SANC, in which Ca 2+ -clock function is preserved, but M clock function is abolished, and therefore AP's cannot occur and do not influence LCR periodicity, clearly demonstrate that: in a fixed, physiologic, free [Ca 2+ ], LCR occurrences are random when phosphorylation level of phospholamban is low; and that LCR periodicity emerges as Ca 2+ cycling protein phosphorylation increases (Sirenko et al., 2014;Sirenko et al., 2013). Prior studies (Lyashkov et al., 2009;Yang et al., 2012) have demonstrated that gradations in the phosphorylation status of phospholamban at Ser 16 across the 5 "steady state" mean APFI's of the present study ( Fig. 9) strikingly resemble gradations of APFI and APFIV across the range of "steady states" observed in the present study (Table 1). At any given time during AP firing, the extent to which clock molecules respond to these activation cues is modulated by their phosphorylation states following a prior AP. bARS or CRs phosphorylate the same proteins that drive the same intrinsic automaticity in the absence of autonomic receptor activation, impacting the memory of synchronization effected by the prior AP.
bARs or CRs affects the synchronization of molecular activation in two ways: directly phosphorylate clock proteins, and indirectly by altering the APFI which alters the cell Ca 2+ activation cues induced by autonomic receptor stimulation.
Our results may provide a cellular basis for a phenomenon demonstrated in a classic study by Nolasco (Nolasco & Dahlen, 1968), i.e., that an AP occurrence, itself, influences the range of APFIs that immediately follow it. The AP, itself, indirectly affects all Ca 2+ -clock functions because it regulates net cell Ca 2+ balance. Functions of M-clock molecules that underlie the generation of an AP indirectly regulate the availability for SR Ca 2+ cycling by modulation of the level of cell Ca 2+ , the SR "oscillatory substrate". Thus, M clock functions also indirectly regulate LCR periods and sizes via their impact on the "steady state" intracellular Ca 2+ level. When the average interval between APs becomes prolonged, a reduction net Ca 2+ influx into: efflux from the cell (Lakatta, 2004) reduces the cytosolic [Ca 2+ ], the rate of Ca 2+ pumping into SR, and the SR Ca 2+ load. These reductions, in turn, prolong the average time from the prior AP occurrence for spontaneous local diastolic ryanodine receptor activation to occur within SANC local micro-domains; the randomness of spontaneous local diastolic ryanodine receptor activation occurring within these micro-domains also increases, broadening the distribution of LCR periods and shifting these to longer times at long AP cycle (Fig. 4).
Thus, the degree of variability in activation of M and Ca 2+ clock molecules that emerges over time following their synchronization by the prior AP is implicated in the cycle length dependence of variability of Ca 2+ and M clock functions measured here (Table 1,  . We interpret the experimentally measured concordant behavior of surface membrane and Ca 2+ regulatory functions during AP cycles within and among the broad range of apparent APFI steady states in the present study to reflect a concordance of degrees of state activation of different molecular species that drive these regulatory functions during an AP cycle. Because membrane and Ca 2+ clocks become coupled in the context of oscillatory electrochemical gradient oscillations during an AP cycle, the extent of self-organized molecular activation within each clock indirectly affects self-organization of molecular activation of the other clock within the coupled system.

Feed-forward self-organization of molecular activation within the coupled oscillator
system. We interpret our results to indicate that both mean APFI and APFIV, the kinetics of AP repolarization (APD90), and the AP-induced Ca 2+ -transient decay time (CaT90), LCR periods, and the kinetics of DD (TNLDD) reflect the extent to which activation states of molecules within the coupled-clock system are synchronized. And further, that mean APFI and APFIV are not only regulated by, but also regulate the degree of synchronization. For example, when an increase in mean APFI reduces net Ca 2+ influx, it indirectly reduces AC-dependent phosphorylation of Ca 2+ cycling proteins, reducing the SR Ca 2+ cycling kinetics and increasing the variability of LCR periods. The prolonged APFI, itself, contributes to the concurrent increase in the AP firing interval variability at long mean APFI. AP characteristics, determined by availability of M clock molecules to respond to a change in membrane potential, both directly and indirectly entrain the Ca 2+ and M clock activation. Thus, mean APFI and APFIV are tightly integrated, and both Ca 2+ and clock molecular activation participate in defining the characteristics of the AP and the timing of the next AP.

Cycle-to-cycle variability of ion channel currents activation.
While cycle-to-cycle variability of LCR characteristics that reflect spontaneous activation of ryanodine receptors can be measured in SANC during spontaneous AP firing, neither ion current amplitudes, nor cycle-tocycle variability can be directly measured in SANC during AP firing. Application of numerical modelling to a time-series of AP clamped experimentally measured APs quantitatively simulated to cycle-to-cycle behavior of ion currents that contribute to APFIV. Numerical simulations of selected ion currents showed that the cycle-to-cycle variability of some current amplitudes informed on the experimentally measured APFIV, while that of other currents did not. For example, under control conditions, cycle-to-cycle variability of numerically stimulated IKr, and ICaL were significantly lower than the cycle-to-cycle variability of the mean APFI and, further, these current amplitudes did not differ significantly between control and CCh (Fig.7D). The moderate cycle-tocycle variabilities of ICaL peak amplitude and very low cycle-to-cycle variability of IKr peak amplitude in control, likely reflect their high importance for governing AP shapes, but relatively less direct importance for variability in diastolic depolarization events that contribute to APFIV.
The numerically simulated INCX decreased (Table 4) when an increase in mean APFI was induced by CCh, in line with the coupled-clock theory of pacemaker automaticity V. A. Maltsev & Lakatta, 2009). Our new finding was that the cycle-to-cycle variability in INCX amplitude in response to CCh paralleled of the APFI (Fig.7D and Table 4). The increase in diastolic INCX indirectly, but explicitly, informs on the evolution of diastolic ensemble Ca 2+ LCR signal. Cycle-to-cycle variability in INCX reflects variability in clock coupling fidelity or effectiveness, because cycle variability of ensemble LCR signal drives cycle-to-cycle variability of INCX amplitude. Cycle-to-cycle variability of LCR shifts concordantly with the mean cycle length across the entire physiological range of APFI (Fig. 3~5, Table 1). Thus, cycle-to-cycle variability of INCX amplitude, not only shifts between control and CCh, as indicated by Fig.7, but also likely across the entire range of "apparent" AP steady states (Table 1 &  In contrast to INCX, cycle-to-cycle variability of If amplitude substantially exceeded that of the APFI in control and became even greater in response to CCh (Fig.7D). This increase in If amplitude between control and CCh (Table 4) is in line with previous ideas regarding the importance of If as an AP cycle "stabilizer" (Noble, Denyer, Brown, & DiFrancesco, 1992), because during longer cycles in CCh, If has more time to activate and reaches a higher amplitude that tends to shorten the very same cycle. It is important to note that while previous patch clamp studies using standard (reductionist) voltage-clamp protocols demonstrated inhibition of If (DiFrancesco, Ducouret, & Robinson, 1989;DiFrancesco & Tromba, 1988) and ICaL (Lyashkov et al., 2009) in the presence of muscarinic receptor activation, our numerical model simulation results were obtained under AP clamp using experimentally measured AP shapes and rates. Thus, our results demonstrate physiological regulation (in the contest of a system's biology approach) of ion currents, i.e. as AP rate physiologically decreases, If has more time to activate (as we mentioned) and ICaL has more time to undergo removal of inactivation in order to reactivate, explaining why ICaL remained unchanged but funny current even increased in our simulations, despite both currents appear to be inhibited when measured by the standard voltage clamp protocols within the reductionist's biophysical approach.
We also found a low variability of IKACh amplitude, similar to that of ICaL and IKr (Fig.7 & Table 4). Low variations of IKACh can be explained by the fact that it is mainly defined by AP overshoot amplitude that remains rather stable. Thus, our analyses demonstrate that the ion current parameters related to the AP phase of the cycle, i.e. peak amplitudes of IKACh, ICaL, and IKr, remain relatively stable (CV < 3.5%), even in the presence of cholinergic receptor activations, whereas parameters related to diastolic phase, INCX (-45 mV) and If, undergo substantial cycle-to-cycle variations concordant with the variation of the APFI. Within the contest of coupled-clock theory , this result can be interpreted to indicate that parameters that exhibit concordant variations with APFI are linked to the time domain, i.e. the diastolic depolarization, whereas parameters that remain relatively stable are linked to the AP phase which must be relatively stable to provide a robust reset to the system clocks.

Long-range beat-to-beat concordance of functions regulate-clock coupling. Long-range
beat-to-beat concordance of functions that regulate clock coupling during each AP cycle in individual, isolated SANC, as measured in the present study, evidenced by principal components analyses (Table 3) and Power Law behavior (Fig. 6), inform on states of order that occur thorough out nature (Bak, 1996) .
Self-similar or fractal-like beating rate variability among cells in vitro has been previously identified in a number of studies in culture but only when cells were confluent, and 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 & DeHaan, 1979;Jongsma, Tsjernina, & de Bruijne, 1983;Kucera et al., 2000). We interpret the self-similarity (power law behavior) among the means of functions measured during AP cycles in different single, isolated SANC operating over a wide range of different apparent "steady states" (Fig. 6) to reflect concordant gradations of self-organized order in clock molecular activation among different cells across the entire physiologic steady states AP firing: the shorter the "steady state" mean APFI, the higher the degree of order (self-organized activation) of clock molecular activation, the more ordered the aggregate of these kinetic functions about that mean (reduced APFIV) (Fig. 5).
The break point in the bi-fractal-like power law relationship (Fig. 6), observed in the frequency domain power law plots, occurs near the transition from the linear to the more nonlinear part of the mean APFI vs. APFIV plot, the time domain (Fig.5). Importantly, the prominent break point in the bi-fractal power law function (Fig. 6) and in distributions of mean APFI and APFIV in the 5 steady states is approximately recapitulated in the relative phosphorylation states of Ca 2+ clock molecules across the range of these steady states (Fig. 9).

Implication of APFIV in single cell to Beat-to-beat variability of AP firing of intact SAN
tissue. Because of added complexity of additional extrinsic modulatory influences on the behavior of SANC imbedded within SAN tissue, concordant graded degrees of intrinsic molecular order shared by individual, isolated SANC as characterized in the present study, cannot directly translate into the mean AP firing rate or its variability of SAN tissue, or into the HR or HR variability in vivo (Bychkov et al., 2020;Jalife, 1984;Michaels, Matyas, & Jalife, 1987;Yaniv, Ahmet, et al., 2014;Yaniv et al., 2011). Electric field potential and phasic membrane, variable autonomic receptor stimulation of SANC embedded within SAN tissue must also potently modulate the mean and beat-to-beat variability of the initiation of intra-SA nodal conduction and exit of impulses that emerge from SAN tissue in vivo (Rocchetti et al., 2000;Yaniv, Ahmet, et al., 2014;Zaza & Lombardi, 2001). But notwithstanding, their modulatory effects in intact SAN tissue, neither the complexity of electrical and mechanical interactions among cells comprising the pacemaker cell syncytium within intact SAN tissue, nor differential modulation of these functions by local difference in autonomic input to SANC tissue can explain the concordant self-similar behavior of clock functions that is shared among a large number of different individual SANC in isolation that is demonstrated by the present results. These intercellular and autonomic factors that regulate intrinsic APFI and APFIV in isolated SANC, however, likely modulate the concordance of selfordered shared, intrinsic molecular activation within individual SANC embedded within SAN tissue. Thus, beat-to-beat gradations in feed forward self-organization of molecular functions intrinsic to SANC, must impact on the beat-to-beat variability of initiation, intra-nodal conduction and emergence of electrical impulses from the SAN (Bychkov et al., 2020;Yaniv et al., 2013a). It follows, therefore, that neither HR nor HR variability in vivo, either directly or specifically, per se, precisely inform on autonomic input to the SANC (Ardell et al., 2016;Billman, 2009). Figure 1. APFI and CaTFI are highly correlated in cells in which both Ca 2+ and membrane potential were simultaneously measured (9 cells, 69 beats).

Figure 2. A-C representative examples of AP recordings measured via perforated patch clamp
(top row) prior to and during bARs or CRs in the same cell; D-F, the density estimates of AP parameters of the self-control data (A-C). The density estimates are nonparametric kernel estimates of probability density functions (Silverman, 1986), and are scaled so that the total area within each curve is unity; G-I, Poincaré Plots demonstrating concordant beat-to-beat variability of AP parameters in control and concordant shifts within the same cell in response to bARs or CRs (the same data from A-C). Note the similar trends between in both fresh and cultured cells.  Table 1. The density estimates are nonparametric 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 5 experimental groups in control and shift concordantly in response to autonomic receptor stimulation in all cells measured.  Table 1. The density estimates are nonparametric 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 5 experimental groups in control and shift concordantly in response to autonomic receptor stimulation in all cells measured.  Table 1prior to and during autonomic receptor stimulation. B.
Surface membrane ion channels, electrogenic ion exchange proteins, e.g., NCX, and ion pumps, e.g., Na-K ATPase are ion current oscillators that regulates TNLDD  and APD90. The NCX and Na-K ATPase are voltage-and transmembrane gradient-dependent, and actively and exquisitely responsive to effect changes in Na + -Ca 2+ electro-chemical gradient oscillations that underline each AP cycle in isolated SANC . Molecules within the bi-directionally coupled oscillator system either directly or indirectly regulate both surface membrane potential and intracellular Ca 2+ . Variable degrees of self-organized coherence or synchronization among these molecular functions impact on the fidelity to which the Ca 2+ oscillators within the system and coupled to the system current oscillators (Y. Yaniv, A.E. Lyashkov, & E.G. Lakatta, 2013b;Yaniv, Lyashkov, et al., 2014).
The SR oscillates intracellular Ca 2+ , that includes the LCR (indicated as the white arrows on Ca 2+ line-scan image) period and CaT90: SR operates as a Ca 2+ capacitor, its Ca 2+ charge is regulated by Ca 2+ ATPase (Serca) that pumps Ca 2+ into the SR lumen, and by ryanodine receptors (RyR) that dissipate the Ca 2+ charge, via releasing Ca 2+ beneath the cell surface membrane. An interaction of phospholamban with Serca modulates the maximum speed of Ca 2+ pumping into SR.
Ion pumps that operate within each clock are energy dependent. The LCR period is an index of clock coupling that becomes manifest when a SANC fires an AP. Figure 9. A Comparison of gradations in the phosphorylation level (A) indexed by phosphorylated phospholamban (PLB) at Ser 16 normalized to total phospholamban of cells in different steady states (redrawn from published previously data (Lyashkov et al., 2009;Yang et al., 2012)) and the means (B) and SDs (C) of all 6 parameters in all 5 "steady states" listed in Table 1. Phosphorylation data and functional data are normalized to value of f-SANC in control, as indicated by the dashed blue line in figure panels (ANOVA F statistic for group effect =23.8, p<0.0001).  recordings were matched by Z-scores to APFI measured by perforated patch clamp. : the cells are the same in which AP were experimentally measured.