A Universal Scaling Relation for Defining Power Spectral Bands in Mammalian Heart Rate Variability Analysis

Background: Power spectral density (PSD) analysis of the heartbeat intervals in the three main frequency bands [very low frequency (VLF), low frequency (LF), and high frequency (HF)] provides a quantitative non-invasive tool for assessing the function of the cardiovascular control system. In humans, these frequency bands were standardized following years of empirical evidence. However, no quantitative approach has justified the frequency cutoffs of these bands and how they might be adapted to other mammals. Defining mammal-specific frequency bands is necessary if the PSD analysis of the HR is to be used as a proxy for measuring the autonomic nervous system activity in animal models. Methods: We first describe the distribution of prominent frequency peaks found in the normalized PSD of mammalian data using a Gaussian mixture model while assuming three components corresponding to the traditional VLF, LF and HF bands. We trained the algorithm on a database of human electrocardiogram recordings (n = 18) and validated it on databases of dogs (n = 17) and mice (n = 8). Finally, we tested it to predict the bands for rabbits (n = 4) for the first time. Results: Double-logarithmic analysis demonstrates a scaling law between the GMM-identified cutoff frequencies and the typical heart rate (HRm): fVLF-LF = 0.0037⋅HRm0.58, fLF-HF = 0.0017⋅HRm1.01 and fHFup = 0.0128⋅HRm0.86. We found that the band cutoff frequencies and Gaussian mean scale with a power law of 1/4 or 1/8 of the typical body mass (BMm), thus revealing allometric power laws. Conclusion: Our automated data-driven approach allowed us to define the frequency bands in PSD analysis of beat-to-beat time series from different mammals. The scaling law between the band frequency cutoffs and the HRm can be used to approximate the PSD bands in other mammals.


INTRODUCTION
The HR is controlled by dynamic and chaotic processes, and it oscillates at different periods over continuously shifting time scales (Goldberger et al., 1990). Therefore, even under resting conditions, mammalian electrocardiographic (ECG) recordings exhibit complex beat-to-beat variations (Billman et al., 2015b). Analysis of the beat-to-beat variation of the HR, a field of research known as HRV analysis, can reveal information about the underlying physiological processes that control its dynamics. For example, a loss of complexity in HRV has been documented in several cardiovascular diseases and has been correlated with an increase in morbidity and mortality (Hillebrand et al., 2013), while abnormal beating patterns can be characteristic of arrhythmias such as atrial fibrillation (Behar J.A. et al., 2013).
Power spectral density (PSD) analysis of the HR fluctuation provides a quantitative non-invasive means for HRV and for assessing the function of the cardiovascular control system (Akselrod et al., 1981). In healthy humans, the frequency of the PSD performed on a 5-min long recording is traditionally divided into three main bands (Malik et al., 1996): the high frequency (HF) band, 0.15-0.4 Hz, where the dominant HF peak can typically be found around 0.25-0.3 Hz; the low frequency (LF) band, 0.04-0.15 Hz, where the dominant peak can typically be found around 0.1 Hz; and the very low frequency (VLF) band, 0.0033-0.04 Hz. The HF band corresponds to rhythms with periods between 2.5 and 7 s (McCraty et al., 2001) and is known to reflect parasympathetic (vagal) activity and respiratory sinus arrhythmia (Akselrod et al., 1981). The LF band corresponds to rhythm modulations with periods between 7 and 25 s and is believed to mainly reflect baroreflex activity while at rest (Malliani, 1995). The LF was also often assumed to have a dominant sympathetic component, but this assumption has been strongly challenged (Billman, 2013b). The VLF corresponds to rhythms with periods between 25 and 300 s. VLF power has been strongly associated with all-cause mortality in cohorts with cardiac failure or multiple organ dysfunction syndrome (Tsuji et al., 1996;Hadase et al., 2004;Schmidt et al., 2005). The energy contained in this band has been suggested to be intrinsically generated by the heart (Kember et al., 2000(Kember et al., , 2001. These frequency bands for humans were standardized by the Task Force of the European Society of Cardiology (Malik et al., 1996) based on the review of years of empirical evidence from various studies (Akselrod et al., 1981;Malik et al., 1996), which applied frequency analysis of the HRV and observed regions of interest within the PSD.
Mammals are commonly used for cardiovascular research. Dogs, rabbits and mice have been of particular interest: dogs are physiologically close to humans and thus a reliable experimental model to study cardiac diseases (Hasenfuss, 1998). The rabbit is the smallest mammal with Ca 2+ dynamics similar to humans (Bers, 2002;Terentyev et al., 2014;Morrissey et al., 2017). With the recent advances in genome manipulation technologies, there has been increased interest in using animals with mutations designed to overexpress or knock out genes implicated in human cardiovascular diseases (Thireau et al., 2008). Mice have been of particular interest in that regard (Tzimas et al., 2017;Hook et al., 2018). Because of the physiological differences across mammals, a quantitative approach is necessary for adapting the frequency bands to different mammalian electrophysiological data. Defining mammal-specific frequency bands is necessary if the PSD analysis of the HR time series is to be used as a proxy for measuring the autonomic nervous system (ANS) activity in animal models.
One parameter that can be used to adapt frequency band cutoffs to different mammals is the HR, and one possible mathematical formulation of such scaling law is the power law. The existence of such a law has not been studied in the context of adapting PSD parameters in HRV analysis, specifically between PSD band cutoff frequencies from different mammals and their respective HR m . Identifying such a power law would enable approximating the PSD frequency bands to mammals other than the ones analyzed in this work.
It has been shown that an allometric law exists between the body mass (BM m ) and various biological processes, such as metabolic rate, life span, HR m , and the ECG PR interval (West et al., 1997(West et al., , 1999Noujaim et al., 2004). This scaling is described by the equation Y = Y 0 · BM b , where Y is the biological process and b is the scaling exponent. Whether an allometric scaling law exists between the PSD band cutoff frequencies and the BM m has not been reported. A secondary aim of this research is to explore whether such an allometric law exists.
This work aims to: (1) create an automated data-driven approach to identify appropriate frequency bands for power spectral analysis of the beat-to-beat time series obtained from different mammalian ECG data; (2) use the estimated bands to research a universal power law between PSD band cutoff frequencies and the typical heart rate (HR m ) of different mammals. This power law could be used to approximate the frequency bands for any other mammal; (3) explore whether a power allometric law exists between the BM m and the characteristic PSD bands.
The public access MIT-BIH Normal Sinus Rhythm (MIT-NSR) database was used for human data (Goldberger et al., 2000). This database was chosen because it contains long ECG recordings of humans having no known cardiovascular condition. Thus, it can be used to define the baseline bands for human HRV analysis. The dog experiments (Billman et al., 2015a) were approved by the Ohio State University Institutional Animal Care and Use Committee and conformed to the Guide for the Care and Use of Laboratory Animals (revised 1996) published by the National Academies Press (Washington, DC, United States). The mouse data (Yaniv et al., 2016) were obtained in accordance with the Guide for the Care and Use of Laboratory Animals published by the National Institutes of Health (NIH Publication no. 85-23, revised 1996). Experimental protocols were approved by the Animal Care and Use Committee of the National Institutes of Health (protocol #441LCS2013). The rabbit data (Brunner et al., 2008;Odening et al., 2012) were recorded in accordance with the local guidelines of the institutions and only after approval by the Institutional Animal Care and Use Committee of Pennsylvania State University College of Medicine and the Milton S. Hershey Medical Center, Hershey, PA, United States, and the Institutional Animal Care and Use Committee of Rhode Island Hospital, Providence, RI, United States, in accordance with the NIH Guide for the Care and Use of Laboratory Animals (NIH Publication no. 85-23, revised 1996). All animal data were obtained from healthy, free-moving animals.
Human data were recorded at 128 Hz, dog data at 500 Hz, and rabbit and mouse data at 1 kHz. All the mammals were conscious, had no known cardiac condition, and no drugs were administered previous to the ECG recordings. For the animal databases, we performed peak detection (Behar et al., 2014) to identify the R-peak locations and we computed the RR interval (defined as the time variation between consecutive R-peaks) time series. We excluded from the dataset large segments of noise. We defined noise segments as the inability of a human annotator to clearly identify R-peaks on the raw ECG trace . Only transient sections of noise were left in the database. In addition, the R-peaks were manually corrected when necessary in order to ensure the reliability of the RR time series. Table 1 summarizes the mammalian databases used in this study. Supplementary Tables S1-S3 provide more details on the individual recordings for each mammal.

Processing the RR Intervals
Before performing PSD analysis, the RR intervals must be preprocessed to ensure that only beats from normal sinus node depolarizations are used. In order to filter out non-sinus beats, we used a moving average filter. We chose the window size to be identical to the values used by the PhysioNet HRV toolkit (Mietus and Goldberger, 2017): a moving average filter with a window size of 21 samples and RR intervals exceeding 20% (for humans, rabbits and mice) or 30% (for dogs) of the window's average were removed. The 20-30% thresholds were empirically Data from humans (Goldberger et al., 2000), dogs (Billman, 2013a), rabbits (Brunner et al., 2008;Odening et al., 2012), and mice (Yaniv et al., 2016) were used.
chosen. In humans, a threshold of 20% is usually used (Mietus and Goldberger, 2017) and it was suitable for the rabbit and mouse data. However, with the dog dataset the 20% threshold was too stringent and removed a high number of normal sinus beats. Thus we decided to select a threshold of 30%, which was more suitable for this dataset. The NN time series is defined as the preprocessed RR time series using the moving average filter.

Power Spectral Density Analysis
The upper bound of the HF band (f HFup ) defines the minimal resampling frequency of the NN time series (i.e., f s ≥ 2 · f HFup , following Shannon's criterion). In the literature, f HFup was set to 0.4 Hz in humans (Malik et al., 1996), but there is no clear rationale for this cutoff. To look for this cutoff we computed the PSD up to the maximal frequency meaningful to resolve (f max ), which was set as the frequency corresponding to the 'characteristic shortest RR interval' found in normal sinus rhythm ECG recordings. We defined the 'characteristic shortest RR interval' as the lower bound of the interval containing 95% of the RR intervals from a large population. We used the histograms of the RR intervals computed on our databases to obtain these values. Given f max , we estimated the PSD on the interval [0-f max ] Hz using the Welch periodogram (Welch, 1967). Because Welch PSD estimation requires uniformly sample data, the RR interval time series was resampled at 2.2 times f max . The resampling frequency was chosen such that it satisfies Shannon's criterion.
To perform PSD analysis, a window (i.e., sub-segment) of the NN time series must be selected. This window must be long enough to resolve the VLF band and short enough to assume data stationarity (Malik et al., 1996). The traditional window for humans is 5 min long (Malik et al., 1996). Because no alternative window length has been standardized for dogs and rabbits, we used the same window size as for humans (i.e., 5 min). We used a 3 min window for mice, as suggested in Thireau et al. (2008). We denote as 'VLF' the frequency band starting from 0.0033 Hz for humans, dogs and rabbits [corresponding to the lowest frequency we can resolve using a 5 min window size 1/(60 * 5)], and 0.0056 Hz for mice [corresponding to the lowest frequency we can resolve using a 3 min window size 1/(60 * 3)], up to the upper bound of the VLF band. In summary, non-overlapping 5 min windows were used for human, dog and rabbit data and non-overlapping 3 min windows were used for the mouse data. PSD was computed for each nonoverlapping window. In the instances where the recording was less than the window size length, the PSD was computed on the total recording length available. This happened in particular for the dog data, where a number of recordings were 4-5 min long.
For each non-overlapping window, we normalized the PSD by the total power in order to allow cross-comparison of mammal types. For each normalized PSD, we detected prominent frequency peaks. Given the location of the detected peaks for each window size, we created a histogram of the prominent peak locations for each mammal type (see example in Supplementary Figure S1). We assumed that the histograms were generated from a mixture of Gaussian distributions; thus, we used a GMM (Reynolds, 2015) to estimate the Gaussian parameters that best describe the underlying distribution. Because of the three traditional power spectral bands (VLF, LF and HF), we used three Gaussians for the GMM (thus assuming three modes). We estimated the GMM and defined the intersection between consecutive Gaussians as the band cutoff frequencies. To set the minimal peak height threshold (defined as the minimal amplitude a peak must have in the normalized PSD to be considered as 'prominent'), we used the human data. Given the selected minimal peak height threshold for humans, we applied the same algorithm to the dog, rabbit and mouse data. We defined the upper bound of the HF (f HFup ) band as three standard deviations away from the mean of the Gaussian describing the HF band. Note that f max defined the upper frequency bound up to which we compute the PSD, whereas f HFup corresponds to the final upper cutoff frequency of the HF band as estimated from the GMM approach. A block diagram summarizing the main steps to identify the frequency bands is shown in Figure 1. We trained the algorithm on a database of human ECG recordings (n = 18) and validated it on databases of dogs (n = 17) and mice (n = 8). Finally, we tested it to predict the bands for rabbits (n = 4) for the first time.

Power-Law
Given the bands identified by the GMM for humans, dogs, rabbits, and mice, we looked for a power-law relationship between the band cutoff frequencies or mean of each Gaussian and HR m . HR m was defined as the median HR evaluated on our database for a given mammal. To search for power laws between band cutoff frequencies and HR m , we used a doublelogarithmic analysis of the band cutoffs for each mammal type against HR m . To search for a power law between the dominant peak in each of the GMM-identified band and HR m , we used a double-logarithmic analysis of the dominant peak location for each 5-or 3-min segment against the mean HR of the segment. We used linear regression to explore whether a linear relationship existed between the variables in the double-logarithmic plot.

Power Allometric Law
To explore this, we used a double-logarithmic analysis of the band cutoffs for each mammal type against the typical BM m of the mammals included in our database. The typical BM m values for the different mammals were taken from Noujaim et al. (2004). We used linear regression to explore whether a linear relationship existed between the variables in the double-logarithmic plot.

Power Spectral Band Definition
We used the human ECG data to train our GMM. The median RR interval for humans was 766 ms and the lower bound of the RR distribution (thus describing the 'fastest' heart rate) was 359 ms (see Table 1). Therefore, f max was set to 2.8 Hz. Based on GMM fitting, the VLF band was identified to be between 0.0033 and 0.046 Hz, the LF band between 0.046 and 0.158 Hz, and the HF band between 0.158 and 0.588 Hz. Figure 2A shows the GMM estimated by overlaying the estimated Gaussians onto the histogram of the prominent peak locations.
We next validated our results on data from dogs and mice. Using the GMM model trained on the human data, the bands were estimated for the dog data; the VLF band was found to be between 0.0033 and 0.067 Hz, the LF band between 0.067 and 0.235 Hz, and the HF band between 0.235 and 0.877 Hz. Figure 2B shows the GMM estimated by overlaying the estimated Gaussians onto the histogram of the prominent peak locations.  (Malik et al., 1996), dog (Billman, 2013a;Billman et al., 2015b), and mouse (Thireau et al., 2008). For dogs, no VLF-LF cutoff could be found in the literature, so only the LF-HF and f HFup vertical lines from the literature are shown. The intersection between two consecutive Gaussians defines the band cutoffs.

TABLE 2 | Gaussian mean (µ) and standard deviations (σ).
Human Using the GMM model trained on the human data, the bands were estimated for the mice data; the VLF band was identified to be between 0.0056 and 0.152 Hz, the LF band between 0.152 and 1.240 Hz, and the HF band between 1.240 and 3.471 Hz. Figure 2C shows the GMM estimated by overlaying the estimated Gaussians onto the histogram of the prominent peak locations.
Using the GMM approach trained on the human data, the bands were estimated for the rabbit data; the VLF band was identified to be between 0.0033 and 0.088 Hz, the LF band between 0.088 and 0.341 Hz, and the HF band between 0.341 and 1.155 Hz. Figure 2D shows the GMM estimated by overlaying the estimated Gaussians onto the histogram of the prominent peak locations.
In addition to the frequency band cutoffs, we also report the Gaussian means and standard deviations for each band of each mammal ( Table 2). Finally, we computed the power ratio: the relative power of each band over the total power for each mammal ( Table 3). Our data show that the HF power is relatively higher in larger mammals (humans and dogs) than in smaller mammals (rabbits and mice).  Figure 3A, we plotted ln(f VLF−HF ) against ln(HR m ).The cutoff frequency scales with the HR m following the power law: f VLF−LF = 0.0037 · HR 0.58 m . In Figure 3B, we plotted ln(f LF−HF ) versus ln(HR m ). Thus, f LF−HF scales with HRm as f LF−HF = 0.0017 · HR 1.01 m . In Figure 3C, we plotted ln(f HF up ) versus ln(HR m ). Thus, f HF up scales with HR m following the power law: f HF up = 0.0128 · HR 0.86 m . Second, we used our GMM approach to determine whether a universal scaling relation exists between HR m and the mean of the Gaussian describing each band ( Table 2). From the linear fit (Figures 3D-F), the following power-law relationships were found: G VLF = 0.0027 · HR 0.53 m , G LF = 0.0077 · HR 0.59 m and G HF = 0.0016 · HR 113 m . Third, we tried to determine whether a universal scaling relation exists between the HR and the dominant PSD peak in each of the frequency bands. In Supplementary Figure S2 we plotted ln(peak band ) for each frequency band (VLF, LF and HF) versus ln(HR). Our data indicate that only the peak HF is correlated with the HR. The dominant peak in the HF band scales with HR as peak HF = 0.0019 · HR 1.09 .
Fourth, we tried to determine whether a universal law for allometric scaling in biology also exists between BM m and the band cutoff frequencies identified by the GMM. We first checked that we could retrieve with our data the known -1/4 allometric law between HR m and BM m as reported in West et al. (1997). Supplementary Figure S3 illustrates the result of the line fit for this analysis. The scaling power found on our data was −0.24, which is very close to −1/4, as expected. From the linear fit (Figures 4A-C)  Finally, we tried to determine whether a universal law for allometric scaling in biology also exists between BM m and the mean of the Gaussians describing each band ( Table 2). In Figure 4D, we plotted G VLF against ln(BM m ) and an allometric law of G VLF = 0.0493 · BM −0.13 m was found. In Figure 4E we plotted ln(G LF ) against ln(BM m ) and an allometric law of ln(G LF ) = 0.1999 · BM −0.14 m was found. In Figure 4F we plotted Frontiers in Physiology | www.frontiersin.org

DISCUSSION
Our first major contribution is the introduction of a new datadriven approach to identify the frequency bands in PSD analysis of heartbeat interval time series from different mammals. To the best of our knowledge, this is the first attempt to provide an automated and data-driven approach for defining the frequency bands in PSD analysis of mammalian HRV. For the human data, our algorithm identified the bands as: Our second major contribution is the application of the method to dog, rabbit and mouse ECG data. Using ECG data from healthy and awake animals (Figure 2), we showed that the frequency cutoffs between the VLF-LF and LF-HF bands could be defined. We validated our model by comparing its results to bands experimentally identified in the literature for dogs and mice: the cutoff between LF and HF in dogs data is similar to the one used in the work of Billman (2013a). Similar cutoffs to ours are also documented for mice: the cutoff between VLF and LF and between LF and HF is similar to other works (Ishii et al., 1996;Just et al., 2000;Joaquim, 2004;Tankersley et al., 2004;Campen, 2005;Fazan, 2005;Adachi et al., 2006;Farah et al., 2006;Baudrie et al., 2007;Thireau et al., 2008). Finally, we tested our model on the rabbit dataset. One study (Goldstein et al., 1995) has looked into defining frequency bands in rabbits. However, this study was performed on anesthetized rabbits, which may require different band definitions than the ones for conscious rabbits, particularly because of the changes in HR that can be caused by the different agents used and the depth of anesthesia (Picker et al., 2001;Mazzeo et al., 2011). Thus, we defined for the first time the HRV frequency bands for this mammal. The prominent peak found in the HF band is characteristic of the respiratory rate (respiratory sinus arrhythmia) (Akselrod et al., 1981). Therefore, we expect the mean of the Gaussian describing the HF band to fall within the respiratory range of the corresponding mammal. The values identified by the GMM (0.397, 0.564, and 2.505 Hz for dogs, rabbits and mice, respectively) fall within the respiratory rate range for these three mammals: 20-40 breaths per minute (brpm) for dogs [0.33-0.67] Hz, 30-60 brpm for rabbits [0.5-1] Hz, and 60-220 brpm for mice [1-3.67] Hz (University Wisconsin-Madison, 2018). Thus the GMM approach successfully retrieved the vagal activity enabled by respiratory effort and which is known to manifest in the HF band.
Our third major contribution is the definition of the upper bound of the HF bands. The cutoff values found in the literature are based on empirical observations. We defined the upper bound of the HF band as being three standard deviations away from the mean of the Gaussian representing the HF band. The upper boundary of the HF band for humans was 0.588 Hz. This is higher than the recommended one (0.4 Hz) (Malik et al., 1996). This higher-than-recommended boundary is likely due to our use of the entire 24 Holter ECG in the MIT NSR database. Indeed, the 24 Holter ECG might include periods of exercise in which the breathing rate might be higher than 24 brpm, i.e., 0.4 Hz. It thus seems reasonable to allow a higher HF upper bound than the standard recommended cutoff for humans. Using the same method, the cutoff frequencies defined for the other mammals were 0.877, 1.155, and 3.471 Hz for dogs, rabbits, and mice, respectively.
Interestingly, our results show that the mouse may serve as a better mammalian model than do the dog or rabbit for studying the effects of drugs, mutations, or cardiac diseases on vagal activity as reflected in the HF band. Indeed, the Gaussian fittings for mice show minimal overlap between the LF and HF bands (Figure 2) compared to other mammals. This means that the vagal activity in mice can be studied without interference from the physiological processes echoed in the low frequency band. However, this is moderated by the larger HF power ratio for humans and dogs (19 and 40.7% respectively) versus rabbits and mice (10.4 and 11% respectively, Table 3). This means that rabbits and mice display relatively less vagal activity than humans and dogs. The effect of respiration is dominant in the HF band. However, during periods of slow respiration, the resulting vagal activity will modulate the HR at frequencies which will cross over into the LF band (Ahmed et al., 1982;Tiller et al., 1996). Thus, for a breathing rate lower than 9.6 brpm, the characteristic vagal frequency peak will fall in the LF band in humans. We interpret this observation as being due to the higher breathing rates of smaller mammals, which leads to a better separation between the respiratory sinus arrhythmia modulation of the heart rate [reflected mainly in the HF band (Supplementary Figure S4)] and other autonomic nervous system effects, partly reflected in the LF bands.
Our fourth major and most important contribution is our discovery of a power-law relationship between band cutoff frequencies, dominant PSD peak, and HR m . Thus, a relationship exists between our GMM-identified cutoff frequencies and HR m . This scaling law can be used to approximate the PSD bands for other mammals not included in our study. In order to define the bands for non-human mammals, other works (Yaniv et al., 2014) scale the frequency bands used in humans linearly with the ratio of human HR m and HR m of another mammal. Figure 3 shows that this approach is incorrect because the relationship between the frequency cutoffs and the HR m is not linear with respect to the ratio of the mean heart rates, as the equations for the log-log subplots highlight. Interestingly, a scaling relation could only be found between the dominant HF peak and HR (Supplementary Figure S2). Because the HF peak represents vagal activity and because the scale power was 1, our results suggest that vagal activity in different mammals is linearly associated with HR. The prominent peak found in the HF band shifts with changes in the respiratory rate (Akselrod et al., 1981), and thus a scaling law was shown between the respiratory rate and the typical HR across mammals. We tested whether the power law that we found can be used to approximate the HRV band cutoff frequencies from other species than the ones included in this work. Aubert et al. (1999) have investigated frequency bands in experimental rat models using direct manipulation of autonomic parameters through pharmacological intervention. They identified the LF and HF bands as: LF [0.19 -0.74] Hz and HF [0.78 -2.5] Hz. The mean HR of the rat population used by Aubert et al. (1999) was 345 bpm. Using this typical HR value and the power law we established between HR m and the PSD cut-off frequencies, we predicted the bands as: VLF [0.0033 -0.110] Hz, LF [0.110 -0.622] Hz and HF [0.622 -1.949] Hz, which is close to the experimental findings of Aubert et al. (1999).
Our fifth major contribution is our discovery that the band cutoff frequencies in mammals, f LF−HF and f HFup , follow an allometric law that scales as the ∼ −1/4 power of the body mass and that G VLF and G LF scale as the ∼ −1/8 power of the body mass. Interestingly, the 1/4 scaling power has been shown to be an essential component in biological scaling (West et al., 1997). In the model of West et al. (1999), processes that rely on hierarchical networks for resource distribution are identified to scale with BM n/4 m (with n ). In practice, many variables in mammalian physiology have been shown to scale with BM m following such a quarter-power law. These include metabolic rate, circulation time, HR, aortic diameter, and respiratory rate (West et al., 1997;Noujaim et al., 2004). Because the frequency bands correspond to the regulation of specific physiological processes (such as baroreflex or vagal activity) on the HR dynamic, it was interesting to find that the band cutoff frequencies (f LF−HF , f HFup ) and the BM m also followed this universal quarter allometric scaling to ensure optimal autonomic activity in regulating the HR dynamics.

Limitations
Although the total length of the recordings was acceptable for the different mammals (Table 1), the sample size (i.e., the number of animals) was limited for mice and rabbits. Moreover, the data were obtained here for specific species of mammals and breeds; however, the scaling laws found can be used to approximate the bands given the typical HR m of another species. For a more precise estimation of the bands in other species/breeds, the same GMM approach can be repeated. In addition, we used data from healthy, free-moving animals. Thus, it might be necessary to adapt the parameters of this analysis for animals during exercise, after drug injection, anesthesia, if genetically modified, or during sleep. Because, the GMM approach is generic, the bands can be easily redefined for other mammalian species or breeds by using the same approach on a representative dataset of animals in these conditions. The ULF band (below 0.003 Hz for humans) has been studied for long term ECG recordings. However, in this work we only consider short segments of no longer than 5 min. As a result, we cannot resolve frequencies below 0.0033 Hz [=1/(60 * 5)]. For that reason, we defined the VLF band as starting at 0.0033 Hz for humans, dogs, and rabbits, and 0.0056 Hz [=1/(60 * 3)] for mice. Finally, a better-adapted window size could improve the accuracy of the bands when using the GMM approach, particularly in the case of dogs and rabbits, where this problem has not been studied. However, when the GMM approach was tested with a 3 min or 5 min window on the mouse data, the cutoff frequencies between the bands were in the same range: f VLF−LF = 0.18 and f LF−HF = 1.33 with a 5 min window versus f VLF−LF = 0.15 and f LF−HF = 1.24 with a 3 min window.

AUTHOR CONTRIBUTIONS
JB, AR, and YY conceived and designed the research. GK, KM, GB, and OS recorded and formatted the data. JB and AR performed data analysis. JB and YY drafted the manuscript. JB, YY, AR, GB, OS, GK, and KM critically revised the manuscript. All the listed authors qualify for authorship and approved the final version of the manuscript. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.