Beating Rate Variability of Isolated Mammal Sinoatrial Node Tissue: Insight Into Its Contribution to Heart Rate Variability

Background Because of the complexity of the interaction between the internal pacemaker mechanisms, cell interconnected signals, and interaction with other body systems, study of the role of individual systems must be performed under in vivo and in situ conditions. The in situ approach is valuable when exploring the mechanisms that govern the beating rate and rhythm of the sinoatrial node (SAN), the heart’s primary pacemaker. SAN beating rate changes on a beat-to-beat basis. However, to date, there are no standard methods and tools for beating rate variability (BRV) analysis from electrograms (EGMs) collected from different mammals, and there is no centralized public database with such recordings. Methods We used EGM recordings obtained from control SAN tissues of rabbits (n = 9) and mice (n = 30) and from mouse SAN tissues (n = 6) that were exposed to drug intervention. The data were harnessed to develop a beat detector to derive the beat-to-beat interval time series from EGM recordings. We adapted BRV measures from heart rate variability and reported their range for rabbit and mouse. Results The beat detector algorithm performed with 99% accuracy, sensitivity, and positive predictive value on the test (mouse) and validation (rabbit and mouse) sets. Differences in the frequency band cutoff were found between BRV of SAN tissue vs. heart rate variability (HRV) of in vivo recordings. A significant reduction in power spectrum density existed in the high frequency band, and a relative increase was seen in the low and very low frequency bands. In isolated SAN, the larger animal had a slower beating rate but with lower BRV, which contrasted the phenomena reported for in vivo analysis. Thus, the non-linear inverse relationship between the average HR and HRV is not maintained under in situ conditions. The beat detector, BRV measures, and databases were contributed to the open-source PhysioZoo software (available at: https://physiozoo.com/). Conclusion Our approach will enable standardization and reproducibility of BRV analysis in mammals. Different trends were found between beating rate and BRV or HRV in isolated SAN tissue vs. recordings collected under in vivo conditions, respectively, implying a complex interaction between the SAN and the autonomic nervous system in determining HRV in vivo.


INTRODUCTION
The normal heart beat dynamics involves orchestration of shortand long-scale periodic signals. These signals are generated by opening and closing of membranal channels (Adair, 2003) in heart pacemaker cells, interaction between pacemaker cells (Michaels et al., 1986), and the pacemaker cell interaction with other body systems (Yang and Xu-Friedman, 2013). To understand the role and relative contribution of each signal, experiments must be performed under in vivo, in situ, and in vitro conditions. When exploring the function of internal pacemaker mechanisms (see for example Yaniv et al., 2015;Behar et al., 2016), the in vitro conditions of isolated pacemaker cells is the optimal experimental model. However, when exploring the interconnected pacemaker cell mechanisms, the in situ environment of isolated sinoatrial node (SAN) tissue isolating it from all environmental effects (hormonal or nervous system) is the ideal model. While ECG recordings (i.e., in vivo) in a variety of mammals and electrical recordings of single pacemaker cells (in vitro) are routinely performed in many labs, electrical data from isolated pacemaker tissue are limited.
The heart rate variability (HRV, refers to variability measured under in vivo conditions) has been suggested as a powerful tool to explore system function (Burg et al., 1993;Bergfeldt and Haga, 2003;Rosenberg et al., 2020). HRV has been quantified in vivo (Goldberger et al., 2000;Behar et al., 2018b) and the beating rate variability (BRV, refers to variability measured under in situ or in vitro conditions) has been quantified in single pacemaker cells (Zaza and Lombardi, 2001;Yaniv et al., 2011). However, although the beating rate of the SAN changes on a beat-to-beat basis, BRV has not been extensively explored in isolated SAN tissue. A number of limitations hinder such research: (i) The electrogram (EGM) is used to measure electrical signals recorded on the isolated tissue surface and reflects the inner currents in this tissue. However, EGM signals differ in beat morphology and rate from in vivo ECG signals even if both are from the same mammal (see Figure 1). Therefore, different analysis tools are required to determine the beating rate from SAN-isolated Abbreviations: ANS, autonomic nervous system; BI, beating interval; BR, beating rate; BRV, beating rate variability (analysis in situ); DFA, detrended fluctuations analysis; ECG, electrocardiogram; EGM, electrogram; FN, false negative; FNR, false negative rate; FP, false positive; FPR, false positive rate; GMM, Gaussian mixture model; HF, high frequency band in the PSD; HRV, heart rate variability (analysis in vivo); IBMX, 3-Isobutyl-1-methylxanthine, a phosphodiesterase inhibitor; IQR, interquartile range, the range between the 25th and 75th percentiles of a dataset; LF, low frequency band in the PSD; MSE, multiscale entropy; PSD, power spectrum density function; SAN, sinoatrial node; TP, true positive; TPR, true positive rate; VLF, very low frequency band in the PSD. tissue EGM than those used to determine beating rates from whole-body ECG recordings. To date, there is no database of mammalian EGM recordings available for the development of such a tool, and there are no standardized, state-of-the-art, partially or fully automated tools to analyze such recordings.
(ii) Assuming that the first limitation is overcome, the beating rate of the tissue can be calculated from the EGM signals. However, there is no standard method to derive BRV from HRV, and there are no publicly available programs to analyze pacemaker tissue BRV. (iii) Isolating pacemaker tissue from healthy human patients is rare; consequently, other mammals are commonly used for cardiovascular research, with rabbits and mice being the most common mammal species used for such research. The rabbit is the smallest mammal with intracellular Ca 2+ dynamics similar to humans (Bers, 2002;Terentyev et al., 2014;Morrissey et al., 2017). On the other hand, mouse models are commonly used for overexpression or knockout of genes implicated in human cardiovascular diseases (Thireau et al., 2008;Tzimas et al., 2017;Hook et al., 2018). Furthermore, mice are practical as aging models due to their short lifespan Yaniv et al., 2016). However, tissues from different mammals differ in their beating rate, and thus, BRV parameters must be adjusted for different mammals (Behar et al., 2018b).
We aim here to overcome these three limitations and design an open-source program to analyze mammal BRV derived from pacemaker tissue EGM recordings. The new tool will be applied to (i) test the effect of drugs on BRV indices, (ii) compare the BRV indices of the different mammals, and (iii) compare BRV indices to their corresponding in vivo indices. This analysis will enhance our understanding of the contribution of pacemaker mechanisms to HRV in vivo.

Databases
EGM data from rabbits (n = 9)  and mice (n = 30) in basal state as well as data from mouse SAN tissues that were exposed to phosphodiesterase inhibition (using 3isobutyl-1-methylxanthine; IBMX) (n = 6) were used . All animal training and validation data used in the present paper were obtained from published studies for which the respective animal protocols and experimental procedures were approved by the local research committee. Rabbit and mouse SAN were fixed in a heated bath (36 ± 0.5 • C) and superfused with Tyrode's solution (see Materials in Yaniv et al., 2014) at a rate of   4 ml/min. An insulated Teflon-coated platinum electrode with a 0.25 (rabbit)-or 0.15 (mouse)-mm diameter tip was placed at the center of the SAN to record extracellular signals using a Neurolog system NL900D (Digitimer, Hertforsdire, United Kingdom), which were recorded at 10 kHz.

Manual Beat Detection
Because no state-of-the-art beat detector is publicly available for mammal SAN EGM data to test our suggested algorithm, beats were manually annotated. The Matlab's "findpeaks.m" algorithm was used for initial peak detection. Then, a single trained annotator reviewed all the recordings and corrected the inaccurate annotations, i.e., false positive and false negative. These reference annotations were then used to evaluate our beat detector and to compute BRV measures. The manual annotations of the training sets were used to calculate the refractory period [minimal beating interval (BI)] and average BI of the SAN from each mammal. Beat Detection Algorithms Figure 2 summarizes the steps used for beat detection in the EGM record, and Figure 3 shows a representative example of the analysis step on one representative signal. In general, each EGM record went through (A) a prefiltering process to clear the data from environmental noise. A notch filter was used, which automatically identifies and reduces effects of the local electricity grid (e.g., 50/60 Hz noise) in the recording. (B) The signal upstream and downstream sign was determined based on the average frequency calculated from the power spectrum of the upward and downward parts of the filtered signal. In order to get more accurate BRV results, the side with the thinner peaks was chosen, reflected in higher average frequency. If the downward direction is preferred, the signal is reversed for the next steps. (C) Naive peak detection: Every 20 s of filtered signal was processed through Matlab's "findpeaks.m" algorithm. The peaks were defined as any point whose distance from a prior beat is longer than the refractory period and were higher than the neighboring points, however wider than 5 data points and not wider than twice the refractory period (width measured at half the height of the peak) and more prominent than a dataderived threshold. Peak threshold = (100 -Q) th percentile of the segment-Q th percentile of the segment (Q = 10 for rabbit and 5 for mouse). (D) After the prominence of all of the naively annotated beats was calculated, those with peaks that were less prominent than 0.7 times the median prominence of beats were eliminated. Finally, the instantaneous BI time series was calculated.

Beating Rate Variability Measures
Prefiltering Before BRV can be calculated, two steps must be taken. First, to assure signal stationarity, a window length of 3 min for mice and 5 min for rabbits was used (Behar et al., 2018b). Second, rangebased filtering was used. A certain constant range was defined, and every beat with a beating interval out of that range was excluded; for this purpose, every BI shorter than the refractory period or longer than three times the average BI of the mammal was discarded (Behar et al., 2018b). The resulting BI time series of the preprocessed signal was named NN.

Time Domain Measures
The majority of time domain BRV measurements is not average BI dependent, thus do not need any adjustment (see Table 1).
As was pointed out before (Behar et al., 2018b), only pNNxx measures that quantify the percent of NN interval differences greater than xx milliseconds must be adjusted for different mammals. We used two approaches to define the xx: one was related to the respiratory rate in vivo, thus the pNNxx in situ was similar to the value in vivo, and the other was to scale the xx parameter according to the scaling ratio of the BRV AVNN relative to the reported HRV AVNN.

Frequency Domain Measures
The Welch's algorithm (Welch, 1967) was used for power spectrum density (PSD) estimation. We chose this spectral estimation method over auto-regressive model (Carvalho et al., 2003;Tarvainen et al., 2006), which is a less frequently used PSD estimate, and over the Lomb method (Lomb, 1976) because of the risk of aliasing (Behar et al., 2018a). Window lengths of 3 min for mice and 5 min for rabbits were used (Behar et al., 2018b). The PSD is traditionally divided into three main bands (Malik, 1996): the very low frequency (VLF) band, the low frequency (LF) band, and the high frequency (HF) band. To determine the cutoff frequencies between the bands, we used a Gaussian mixture model (GMM) of 2 Gaussians on the histogram of all prominent peaks (see Figure 4), following the approach described in Behar et al. (2018a). To calculate the prominent peaks in each band, we used a simple peak detection algorithm to look for the 16 most prominent peaks on each of the normalized PSDs, with a threshold of 0.01. The minimal frequency was determined as one over the window length in seconds and was used to defined the lower band of VLF. The HF band was set to be between the high cutoff frequency of the LF band and 2 Hz.

Non-linear Domain Measures
Four measurements were used in this group: β coefficient, which corresponds to the slope of the linear interpolation of the spectrum in a log-log scale for frequencies below the upper bound of the VLF band, detrended fluctuations analysis (DFA) measures (Peng et al., 1994), Poincare analysis variation measure, and multiscale entropy (MSE) measures (Costa et al., 2005). Poincare analysis and MSE measures require no adaption. The β coefficient was estimated in the adjusted VLF band. Originally, two DFA coefficients were reported for the slopes before and after 16 beats (cutoff). We evaluated this cutoff for BRV analysis.

User Interface
The PhysioZoo open source software (Behar et al., 2018b) was used to calculate the different BRV measures. All of our analysis tools were implemented in PhysioZoo and are open to the public.

Performance Statistics
The quality of the beat detector was assessed by comparing between the results of the algorithm and the manually annotated beats, by using a tolerance window of 10% of the average BI of the dataset and ignoring a constant difference between the beat locations. True positive (TP) annotations, false positive (FP) annotations, and the number of missed beats [false negative (FN)] were calculated. The following statistical measurements were used to report on the quality of our beat detector: True positive rate (TPR)-the percentage of correctly annotated beats out of all the real beats of the dataset False discovery rate (FDR)-the percentage of falsely annotated beats out of all the beat annotations False negative rate (FNR)-the percentage of missed beats (real and not annotated by the algorithm) out of all the real beats of the dataset

General Statistics
The rank-sum test was used to define the significance level of the differences between in situ vs. in vivo conditions of the same animal and between mouse vs. rabbit SAN tissue.

RESULTS
This section presents the results of the performance of the beat detector on mouse and rabbit SAN tissues recordings collected under basal conditions, the range of BRV measures obtained for mouse and rabbit SAN tissues, the interpretation of BRV results in comparison to HRV, and an example of insights on pacemaker function gained from drug response BRV analysis.

Beat Detector
To validate the ability of the beat detector to perform on unknown EGM records, we divided the mouse data into training, validation, and test datasets. The dataset was randomly divided into a training set (40%), validation set (20%), and test set (40%). The algorithm was developed using the training set and then fine-tuned by evaluating its performance on the validation set. Finally, the generalization performance of the algorithm is reported for the test set.
In the case of the rabbit data, because of the limited number of animals, the data were divided into training (67%) and validation (33%) sets. Table 2 provides the performance statistics of the detector for mouse and rabbit. The beat detector algorithm very accurately detected the beat in the SAN EGM. Figure 5A presents the BI histogram of representative examples of human-curated mouse and rabbit data. As can be seen, the BI scattering of mouse SAN tissue was higher than in the rabbit The Poincare plot ( Figure 5B) was more scattered in mouse than rabbit, in accordance with the results of BRV time domain parameters. Figure 6 presents the log PSD vs. log frequency for all data; a tight linear relationship between the log (PSD) and log frequency in the VLF band was achieved. Thus, we can define β as the slope of this linear regression between the log (PSD) and log frequency in the VLF band, similar to that obtained by HRV analysis. Figure 7 visualizes that the cutoff block size of DFA was similar to the cutoff under in vivo conditions [19.3 ± 7.6 (n = 7) for rabbit and 17.5 ± 9 (n = 11) for mouse]. Table 3 summarizes the major changes in BRV analysis under in situ conditions vs. in vivo. Table 4 summarizes the analysis results of BRV measures of rabbit and mouse. The time domain analysis found a lower variability in rabbit vs. mouse SAN tissue. In both mammals, the HF band had the least information about PSD,   as expected. The PSD information was divided between the VLF and LF bands (ratio of 4:1 in mouse and ratio of 1:2 in rabbit between VLF and LF bands). The non-linear parameters showed lower complexity in mouse than rabbit.

BRV Measures
To further explore the variability and complexity observed in SAN BRV, we plotted the Poincare plot and MSE vs. order, respectively. Figure 8 shows that, in the lower scale, the MSE curve in rabbit was lower than that of mouse. However, at higher scales that represent the system complexity, the trends reverse.

Studying BRV in Response to Conditions That Affect Pacemaker Function
To study how direct changes in pacemaker function affect the BRV, phosphodiesterase activity was inhibited by applying 100 µM 3-isobutyl-1-methylxanthine (IBMX),   12). A change in slope is visually noticeable around 15 < n < 20. The traditional value used for ECG-derived HRV is n = 16 (Peng et al., 1995) and shown with a dashed vertical line.  The band cutoff frequencies, the XX calculation for pNNXX, and minimum and maximum BI were calculated using the entire in situ dataset.
which subsequently leads to an increase in phosphate activity and to an increase in BR. After 5 min of exposure to 100 µM IBMX, the BRV time domain measures are decreased (Table 5 and Figures 9A,B). Additionally, a reduction in the relative PSD was documented in the LF, alongside an increase in the relative PSD in the HF band in response to 100 µM IBMX. There was also a shift in the LF peak ( Figure 9C). The MSE curve showed an increase in the complexity in low-scale factors in response to 100 µM IBMX and similar behavior at high scale ( Figure 9D).

Beat Detector
We present here a novel beat detector that is suitable for EGM data. Both the beat detector and the annotated data used to evaluate its performance are available at PhysioZoo platform (Behar et al., 2018b). To the best of our knowledge, this is the first time that a beat detector and annotated database are published. The beat detector was trained on mouse data and tested on a separate database with recordings from different mice. Because the amount of rabbit data was limited, the beat detector was trained and tested on data that came from the same rabbit. Future testing with recordings from additional rabbits and validation and training with other types of mammal data (see limitation) are expected to improve our beat detector.

BRV Measures
We present here, for the first time, a procedure for calculation of BRV measures in different mammals. The main adaptation from HRV to BRV is the changes in the band cutoff frequencies and the XX for pNNxx. Although the BR was lower under in situ conditions, the cutoff frequency between LF and HF was higher. These results are in contrast to the trend observed under in vivo conditions (Behar et al., 2018a). We also justified that the range of BI filtering should be different between in vivo and in situ conditions due to change in average BR, as expected. Under in vivo conditions, larger animals have a lower BR with increased HRV (Noujaim et al., 2007;Behar et al., 2018a). However, in isolated SAN tissue, a larger animal has lower BR but with lower BRV. Thus, the non-linear inverse relationship between the average HR and HRV is not maintained under in situ conditions. Because both the autonomic nervous system (ANS) and the SAN control the HR and HRV, lower HRV in small mammals is achieved via the balance between the sympathetic and parasympathetic systems. Thus, the ANS and the SAN have opposite effects on the HR vs. HRV relationship. The division between the two ranges in DFA is assumed to be related to the relationship between the beating and breathing rates (Willson et al., 2002). However, we found that the cutoff between the two ranges under in vivo and under in situ conditions was similar. Thus, one may doubt that the breathing rate is the cause of the division, and it may be related to some internal pacemaker mechanism.

BRV in Response to Conditions That Affect Pacemaker Function
We further explored the BRV in response to a drug that increases the beating rate (100 µM IBMX). As was shown before , an increase in mouse BR was associated with a decrease in the BRV time domain. Additional data regarding  how drug intervention affects rabbit BRV, and BRV analysis from other mammals, will provide further insight into the BR vs. BRV relationship in situ. Table 4 summarizes the main differences between in vivo and in situ BRV. Based on the time domain analysis, there was an increase in fragmentation parameters in recordings of isolated tissue as compared to in vivo conditions (Behar et al., 2018b).

BRV vs. HRV Measures
Exploring the frequency domain parameters revealed interesting phenomena. The majority of information was restored in the LF and VLF bands of isolated tissue compared to in vivo conditions. Thus, the SAN is the main contributor to the VLF and LF bands in rabbit and to the VLF band in mice. One of the consequences of such behavior in mice and rabbits can be related to the balance between sympathetic and parasympathetic activity and to the relative magnitude of each of these activities in different mammals. The mice in our research were hosted in 20 • C well below the thermoneutrality temperature and thus showed higher sympathetic activity, and sympathetic activity affects LF more than HF (Axsom et al., 2020). However, acclimation to a thermoneutral environment reversed the balance between sympathetic and parasympathetic activity (Swoap et al., 2008;Axsom et al., 2020), and thus, using SAN tissue from such mice in the future can reverse the information restored in different bands. As expected, there was almost no PSD in the HF band in both rabbit and mouse SAN tissue. Thus, the ANS is the main contributor to that band. Note that respiratory rate frequency is characteristic of the vagal activity. Because vagal activity is expressed in the HF band, it can explain the reduction in that band in the isolated tissue.
Based on the non-linear analysis parameters, the system complexity in the rabbit SAN tissue was higher as compared to in vivo conditions; the opposite trend was observed in mice. Thus, the system complexity is maintained mainly by the SAN in rabbit, and in mice, both SAN and ANS contribute to this complexity.
We and others have shown before that smaller mammals have a higher HR than larger mammals (Noujaim et al., 2007;Behar et al., 2018a). We also showed that an increase in HR was associated with a decrease in sample entropy ( Figure 10A). In contrast to the in vivo conditions, in isolated SAN tissue, the opposite relationship between entropy index and BR was observed ( Figure 10B). Thus, putting the in vivo and SAN tissue data together on one graph is expected to show biphasic behavior for the BR and BRV. One can therefore hypothesize that when the beating rate increases in situ, there is an increase in shorttime entropy quantified by first-order entropy. Thus, the intrinsic properties of SAN lead to an increase in entropy in response to an increase in BR. However, the fact that this relationship reverses in vivo implies that (i) the ANS leads to an opposite relationship between entropy and BR, and (ii) under in vivo conditions, the ANS controls this first-order entropy scale as was shown before (Rosenberg et al., 2020).

LIMITATIONS
The beat detector was trained and tested on SAN from healthy animals, without any drug interventions. Future training and examination of the beat detector on EGM of SAN tissue that was exposed to a drug will be necessary to validate the use of the beat detector. Similarly, we analyzed BRV using recordings of SAN of healthy mammals. Future analysis of BRV of SAN tissue from transgenic animals or animals with cardiac diseases and comparison to the respective in vivo conditions will allow us to learn how changes in SAN function affects HRV in vivo.

CONCLUSION
The approach presented here will enable standardization and reproducibility of BRV analysis in mammalian. Different trends were found between BR and BRV of isolated SAN tissue vs. HRV in vivo conditions, implying a complex interaction between SAN and the ANS in determining HRV in vivo.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in PhysioZoo.com.

ETHICS STATEMENT
The animal training and validation data used in the present paper were obtained from published studies for which the respective animal protocols and experimental procedures were approved by the local research committee.