The Fundamental Connection Between Typical Time-Domain, Frequency-Domain and Nonlinear Heart-Rate-Variability Parameters as Revealed by a Comparative Analysis of Their Heart-Rate Dependence

1. Introduction

Heart rate variability (HRV) measures the variation in interbeat intervals from one heartbeat to the next, reflecting changes in the heart rate (HR) over time. Heart rhythm is primarily governed by the synchronized firing of pacemaker cells in the sinoatrial node of the heart muscle, which initiates the cardiac cycle, and the activity of the pacemaker cells is regulated by autonomic efferent neurons and circulating hormones. The complex actions of the underlying interdependent regulatory systems give rise to the variability of the length of the cardiac cycle over different time scales, supporting the optimal performance of the heart under homeostasis. Therefore, HRV is widely considered an important marker of the autonomic nervous system (ANS)[1].

In addition to sympathetic and parasympathetic nervous system (SNS and PNS) activities, respiration and the baroreceptor reflex are considered to influence HRV on time scales ranging from seconds to minutes. In contrast, other ultradian and circadian rhythms driven by factors like metabolism, physical activity, body temperature and sleep-wake cycles govern HRV over hours to the daily cycle.

At the individual level, HRV can be influenced by several non-modifiable factors, such as age, gender and genetics, as well as modifiable factors, including physiological and pathological elements (e.g., posture, respiration, hormonal influences, inflammation, infections, neurological or mental disorders, cardiovascular conditions), lifestyle factors (e.g., alcohol, tobacco, drug use, exercise, meditation) and environmental factors (e.g., physical stress and temperature), and neuropsychological factors [2,3,4,5].

HRV is widely used in the investigation of the heart, and is considered an essential diagnostic tool in cardiology. Higher HRV is often associated with better fitness and heart adaptability, while reduced HRV is commonly linked to various pathological conditions, including congestive heart failure, diabetic neuropathy, mental disorders and cancer [6,7].

Researchers and clinicians use time-domain, frequency-domain, and non-linear indices to measure HRV, by typically analysing time-resolved electrocardiography (ECG) or pulsometry signals [8,9,10]. Time-domain indices are derived from the time series of RR-intervals, and capture the variability in HRV over the monitoring period. The mean RR (or NN) interval or, equivalently, the average HR, are the simplest time-domain indicators. Other time-domain measures, such as the standard deviation of normal-to-normal RR intervals (SDNN), represent the overall variability (both short- and long-term), while the root mean square of successive differences (RMSSD) is more suited to quantify short-term variability [8].

Frequency-domain analysis involves calculating the absolute or relative power, within high-frequency (HF), low-frequency (LF), and sometimes very-low-frequency (VLF) bands. For HRV in normal subjects, typical frequency band ranges in short-term recordings are 0-0.04 Hz (VLF), 0.04-0.15 Hz (LF), and 0.15-0.4 Hz (HF). Frequency-domain HRV measures include peak frequencies (the frequency corresponding to maximum power within each band), as well as the absolute and relative powers for each band, sometimes the normalized powers for LF and HF, the LF/HF power ratio, and the total spectral power [9,10].

Nonlinear indices assess the complexity and unpredictability of the interbeat intervals. Due to the heart's complex regulation system, HRV cannot be fully described using linear methods alone, prompting the use of nonlinear techniques such as the Poincaré plot, sample entropy (SampEn), and detrended fluctuation analysis (DFA). The Poincaré plot visualizes the correlation between consecutive RR intervals by plotting RR_j+1 against RR_j, and the shape of the plot is quantified using an ellipse fitted to the data. The ellipse's width (SD1) and length (SD2) are considered to represent short-term and overall HR variability, respectively [11].

Sample entropy SampEn measures the complexity of HRV, with a computation depending on two parameters: the embedding dimension (m) and the tolerance (r). Default values are m = 2 and r = 0.2xSDNN. Detrended fluctuation analysis (DFA) evaluates correlations in HRV data across different time scales. Short-term (α1) and long-term (α2) fluctuations are measured by the slopes of a log-log plot of correlation coefficient versus segment length, with α1 derived from segments of 4-16 beats and α2 from 16-64 beats [8,9].

HRV analysis is nowadays routinely applied in clinical cardiology and stress-relaxation studies [12,13,14,15]. When choosing HRV measures from the time-domain, frequency-domain, and non-linear methods to track clinical progress or performance, the choice often depends on the investigator's intuition and preference. This choice is further complicated by the well-known but often ignored fact that changes in HR do affect HRV. E.g., a decrease in HR leads to an increase of typical HRV indices, while an increase in HR generally decreases HRV parameters such as SDNN and RMSSD [16,17,18,19]. As physical activity and posture impact both HR and HRV [3,4,5,20,21], protocols developed by the joint European and American task forces [22] recommend stationary RR time-series measurement during HRV monitoring, though this complicates experiments requiring movement. Nearly all parameters that describe instantaneous HRV (i.e., excluding time-averaged measures) are influenced by HR, such as the time-domain, frequency-domain, or nonlinear measures [23,24,25,26,27], and this dependence is generally pivotal, namely, it usually masks the effects of all the other factors on HRV. For instance, if the relationship between the common HRV parameter, RMSSD and the HR is not appropriately considered, then the variability parameter will primarily reflect changes in the heart rhythm, and it will not provide independent physiological information [27].

During the past two decades, some prominent research groups have, indeed, recognized that HRV cannot be used as an accurate diagnostic marker without taking into account its HR-dependence, and therefore serious attempts have been made to describe the general features of the HR-dependence of various HRV parameters (see, amongst others, the publications of and the debates between the Billmann and the Boyett groups [16,16,18,19,23,24,25,28]. However, lacking the exact function that describes the relationship between HRV and HR, there has been no practical way to define a single, ”heart-rate-corrected” HRV parameter that could properly account for the HR-dependence. Some papers even argued that the existence of such a well-defined HRV(HR) function might mean that HRV does not contain extra information to what HR does [24]. Lately, it has been shown that while a well-defined RMSSD(HR) function indeed exists, and shows a remarkable invariability on the level of an individual from the hours’ to the monthly scale (hence it is called the “Master Curve”, MC) [27], momentary deviations from it on the seconds’ to the minutes scale do occur, and reflect changes in the sympatho-vagal balance, due, e.g. to stress or relaxation [29]. Longitudinal analyses also revealed changes of the MC on the decade scale, as well [27], which may be associated with the age-dependent degradation of the cardiac regulatory mechanism. As for the HR-dependence of the spectral and nonlinear HRV parameters, we should highlight the contribution of Platisa and coworkers, who, based on the results of some 20 volunteers, presented nice, low-noise curves describing changes in the VLF, LF and HF components, as well as in the sample entropy and alpha1 and 2 parameters of the DFA analysis, as a function of RR intervals, in a row of publications [26,30,31,32,33]. They also established apparent pairwise connections between these parameters, described with quasi-linear or quasi-quadratic functions, in this representation. More recently, Tsaneva and coworkers also presented notable analyses of HRV data by time-and frequency-domain, as well as nonlinear descriptors [10,34,35,36,37].

All the above works, along with many others, are fundamental contributions towards a deeper understanding of the complex mechanism via which the ANS influences the fluctuations of the HR time series, however, a more complete exploration of this stochastic process is at need.

Correspondingly, the current study aims to address this gap by performing a comparative statistical analysis of ECG recordings from a public database, containing 24-hour-long ECG recordings of 200 healthy individuals, with a focus on the HR-dependence of typical time-domain, frequency-domain, and nonlinear HRV parameters, which reflect the same stochastic process from different viewpoints. in addition, the large number of recordings and the wide age distribution of the volunteers allowed us to make age-cohort studies, as well. In support, the age dependence of the HRV parameters illuminates the physiological grounds of the complex regulatory mechanism of the ANS on cardiac function, since aging is associated with structural and functional alterations in the brainstem and hypothalamus, regions integral to autonomic regulation [38,39],.

In this study, we represent the various HRV parameters as a function of HR, because the curves in this representation are the easiest to interpret in terms of PNS and SNS effects (see, e.g., Figure S2 and S3) on the ultradian time scale of minutes to hours. Age-dependence, on the other hand, reflects long-term changes in the cardiac control of ANS, typically on the scale of decades to years. Investigation of such a combination of dependencies on a large database is a novel approach to the best of our knowledge. Therefore, we propose that our results, clarifying the formal connections between the most typical HRV parameters via analyzing their HR- and age-dependence, will facilitate a more insightful interpretation of HRV data in cardiovascular research and medical practice.

2. Materials and Methods

ECG data were obtained from a public database of the Telemetric- and Holter-ECG Warehouse (“THEW”) at the University of Rochester Medical Center, New York, United States [40]. The 24-h Holter recording data of 202 healthy volunteers (Database Normal, EHOL-03-0202-003, age ranging from 9 to 82 years) were analysed. The time series were filtered for outliers by a MATLAB routine “isoutlier”, with moving median 30 points (MathWorks, 2020).

In this study, we examined time-domain HRV parameters RMSSD and SDNN directly calculated from the time series of RR intervals extracted from the ECG recordings, the frequency-domain parameters LF, HF and very-high-frequency (VHF) power obtained from Fast Fourier Transform (FFT)-spectra of the RR and dRR time series, while among the nonlinear measures the DFA exponent α1 (with m = 2 and r = 0.2xSDNN), and the Poincaré plot were derived from RR time series, according to their definitions (Table 1), using a MATLAB code written by one of us (A.B.). Note that, since here we examined only healthy subjects’ recordings, where the percentage of extrasystoles is negligible, we approximated the NN time series with the corresponding RR time series, filtered for outliers (Figure S1).

For the DFA analysis, the time series was divided into non-overlapping segments of lengths ranging from 10 to 100 data points, in 10 steps. These fixed segment sizes (10, 20, ..., 100) were used to evaluate the fluctuation function across different time scales. The RR interval time series was interpolated at 4 Hz to achieve uniform sampling. All analyses — including frequency-domain measures (FFT-based spectral analysis), nonlinear metrics (e.g., entropy calculation, detrended fluctuation analysis), and time-domain parameters (e.g., RMSSD, SDNN) — were performed on successive overlapping segments. Each segment contained 512 data points, corresponding to a time window of 128 seconds. Consecutive segments were shifted by 50 data points (12.5 seconds), resulting in a considerable overlap between successive windows, which allowed for high temporal resolution in the dynamic analysis.

HRV-related Poincaré plots are graphical representations of the relevant RR time series [11]. The Descartian coordinates of the i-th point are R_i and R_i+1, respectively, where R_i is the i-th RR interval (Figure 1a). The Poincaré plots can then be transformed to the so-called Bland-Altman representation ((R_i + R_i+1)/2, R_i+1 − R_i). Note that the Bland-Altman plot is a difference plot, usually used to compare two different series of measurements of the same thing or process, in order to reveal, eg., possible systematic differences, due to inaccuracies of the different measuring instruments [43]. Here, we used it in an unusual way, inasmuch as we applied it to the same time series, but one of the counterpart is ”shifted” by one data point. and this way, we arrive at a transformation of the Poincaré scatter plot, where the points are now distributed quuasi-symmetrically around the abscissa (Figure 1b and c). From this, after a conversion from the RR- to the HR-scale (Figure 1c), and fitting Gaussians to the dRR distribution at each HR, a high-accuracy RMSSD(HR) plot can be determined, which is called the “Master Curve” (MC) (Figure 1d) [27].

The evaluation methods were applied to each valuable ECG recording, and the relationships between the time-, frequency-domain and nonlinear parameters are shown via HR data of a typical healthy individual (age:21, male), doubled with age-cohort-averaged data of healthy people.

3. Results and Discussion

3.1. RMSSD and the Poincaré plot

Revealing the connection between these time-domain and nonlinear measures requires several-hours-, preferably 24-hours-long recordings. Figure 2a shows a Bland-Altmann-type representation of the conventional Poincaré plot of a typical individual from the THEW data base, where from the MC, representing the HR-dependence of RMSSD, was calculated, as described above (Figure 2b). RMSSD as calculated directly from the RR time series data by a conventional sectional evaluation of the HR time series with a rolling time window (width = 128s, step size = 12.5s) is shown together with the MC in Figure 2c. When selecting the optimal width of the time window, one has to consider the problem stemming from the time dependence of HR. If the analysis is restricted to short time intervals (where HR does not change significantly), it results in a higher uncertainty in HRV, as compared to the MC method (Figure 2c). Although, longer averages reduce the level of fluctuations in the sequential time-window analysis, as well, but they inevitably accompany with information loss on the HR scale. Hence, whenever it can be applied, the MC method is preferred to use for the determination of the RMSSD(HR) dependence. As it was established in Buzás et al [27], the MC is rather stable for an individual on the daily to monthly scales, so it can be used as a good reference for the cardiac state of a person. On the cohort level, it was pointed out, however, that aging and myocardial infarction do influence MC, so a change of MC is considered to reflect remodelling of the heart. In Figure 2d, we show the age-cohort-averaged Master Curves according to a new age classification (the start and end ages of the cohorts are also indicated on the plot), showing a gradual decline of the HRV level, as age progresses [27]. It has also been shown that MC can be decently interpreted by a two-component stochastic model, where the two noise components were tentatively attributed to sympathetic and parasympathetic influence, PNS effects being more dominant at low heart rates (HR< ca.100), while SNS effects above (Figure S3) [27]. Recent publications suggest that MC represents the mean sympatho-vagal balance as a function of HR, and short-term (i.e., minute-scale) deviations from it are indicative of actual fluctuations of the SNS and/or PNS effects [27,29]. A new HRV-indicator based on the actual RMSSD normalized to the Master Curve was successfully used to characterize stress responses and relaxation [29].

3.2. Fourier-components versus RMSSD and SDNN

Fourier-transform of the RR or dRR time series is used for the spectral description of HRV, usually using the FFT algorithm. FFT analysis of the RR(t) data usually show characteristic features in the “very-low-frequency” (0.0033Hz<ω<0.04 Hz), the “low-frequency (0.04Hz<ω<0.15Hz) and in the “high-frequency” (0.15Hz<ω<0.4 Hz) regime, denoted by VLF, LF and HF, respectively. In the FFT spectra of dRR data, however, the VLF components are overdamped, hence, usually only the LF and HF features are considered. The VHF regime (>0.4 Hz) is rarely mentioned and poorly discussed [44].

Unless the so-called ULF band (ω<0.0033Hz) is investigated, FFT is usually performed on HRV(t) (i.e., RR or dRR time series) data within a few-minutes-wide time window (of a width usually between 2 and 5 minutes), which is shifted along the time scale. The RR or dRR time series should be interpolated by equidistant sampling beforehand, and then an average HR value can also be assigned to each interval. Because of the strict mathematical relationship between frequency- and time-domain parameters, if the former show a HR-dependence, the latter should also do so.

Figure 3a shows the averaged FFT map of all the healthy people’s dRR data from the THEW database, as a function of HR. While the conventional LF and HF bands present explicit spectral features at “normal” heart rates (HR up to ca. 120 bpm), a more diffuse but non-negligible feature also emerges in the VHF region at higher heart rates (above ca. 100 bpm, Figure 3a and b). Without entering a deeper discussion about the origin of this neglected broad band, one should establish that it appears in the range of high HRs, where SNS effects dominate, and as such, it is expected to be an important status indicator in neurology, cardiology, aging or sports applications.

Our findings provide an interesting support to the above statement: The 24-hours-long ECG recordings of the THEW database, registered on healthy people of ages between 9 and 82 years, were classified into 5 age cohorts, and the cohort-averaged FFT spectral maps were compared (Figure 4).

It is apparent that as age increases, a large part of the VHF “cloud” is shifted towards lower HR values, and gradually merges into the HF band. Further targeted physiological studies are needed to explore the full diagnostic value of this new spectral feature, but it is clear that this effect points at the increased weight of SNS effects, as age progresses. Our new findings are also in line with independent physiological evidences for the increase of the SNS activity by age in humans, which has been observed in both basal (resting) conditions and in response to stress or physiological challenges. In support, Monahan et al. (2003) demonstrated that aging is associated with increased sympathetic nervous system activity and reduced parasympathetic modulation, which may contribute to the development of age-related hypertension [39]. Esler et al. (2008) observed elevated plasma norepinephrine levels and increased sympathetic nerve firing rates as aging is progressed [45].

Due to the strict connection between the HRV(ω) and HRV(t) data (where ω is the angular frequency, and t is time), a quantitative statement can be established for the HR-dependence of such time-domain HRV parameters as RMSSD or SDNN and that of the Fourier-components of the FFT power spectra. Namely, the following equations hold:

The name of the underlying mathematical relationship is Parseval’s theorem, which expresses the unitarity of the Fourier transform, loosely stating that the “energy” of the signal should be equal by summing power-per-sample across time or spectral power across frequency [46]. To translate it for our case: the HRV amplitude associated for a certain HR should be the same whether it is calculated from the time-domain parameters or the corresponding frequency-domain ones. From this, it follows that RMSSD(HR) and SDNN(HR) can be regained from the superposition of the LF(HR), HF(HR) and VHF(HR) functions. . Using the data of the selected typical individual, Figure 5 shows that it is indeed the case, with good approximation.

3.3. Fourier-components versus nonlinear measures (DFA and SampEn)

DFA analysis and SampEn report rather on the structure of HRV fluctuations than on the size of them. DFA analysis evaluates the degree of correlation between each RR interval and the preceding and following RR intervals across various time scales. The algorithm first divides the integrated time series of 10 intervals into bins of length N. In each bin, a least-squares line is fitted to capture the local trend. The integrated time series is then detrended by subtracting the local trend in each bin. The root-mean-square fluctuation of the integrated and detrended time series is subsequently calculated for all bin lengths, ranging from N = 10 to N = 100 RR intervals. The α coefficient, which quantifies the relationship between fluctuations and bin length on a log-log plot, determines the extent of long-range correlations. An α₁ value of 1.00 indicates perfect correlations across all time scales, while an α value of 0.50 signifies complete lack of correlation with previous and subsequent intervals [8,47].

Figure 6d shows the α₁ (HR) plot of the selected individual, whose RR data were used in the former HRV analyses in this paper. In spite of the relatively large standard deviation, the overall trend is clearly visible: While at the lowest HR values (HR ≈ 60-65 bpm), the DFA coefficient, α1 is close to 0.5, there is a rising tendency up to HR ≈ 120-130 bpm, α1 reaching a maximum value of shortly over 1, above which HR range, it starts to decline to lower values, again. The statistical meaning of this is that while HRV statistics at low HR values is close to that of an uncorrelated white noise, around its zenith it approaches α1≈1, representing a strongly correlated pink (1/f) noise, and at even higher HR values, it returns toward a less correlated fluctuation statistics. Entropy, reporting the degree of disorder, reflects a corresponding anti-correlation with the DFA exponent. The uncorrelated white-noise-like state corresponds to the highest value, and then, after reaching a minimum showing a stronger level of correlation at HR ≈ 130 bpm, it starts to rise again as HR grows (Figure 6a). The corresponding frequency spectra are in line with the above: at HR ≈ 65 bpm, the noise spectrum shows a relatively flat feature (slope ≈ 0), at other HR values it shows a steeper negative slope, reaching its extremum (slope ≈ -2) around 130 bpm, representing an 1/f²–like, correlated noise.

An age-cohort analysis shows the above relations between the DFA-exponent and SampEn even more clearly (Figure 7). The cohort-averaged plots show a nice anti-symmetry, revealing that the decline of α1 at high HRs is real, whose pandant is seen as a rise of SampEn in the same range. Our findings are in concert with the earlier results of Platisa et al., who depicted the DFA-exponent and SampEn as a function of the RR-interval lengths, and found extrema around the same RR value in both cases [30]. They also hypothesized that the abscissa of this “turning point”, H_tp (i.e., the minimum place for SampEn and the corresponding maximum place of α1) might reflect the intrinsic heart rate, H_i (i.e., the HR measurable at full autonomic blockade) [30,33]. This would mean that H_i could be determined non-invasively, without pharmacologically blocking the effects of PNS and SNS on the cardiac cycle. Our age-cohort-averaged data, indeed, show a similar tendency of declination for H_tp, as it has earlier been determined for H_i [48], implying a possible close connection between H_tp and H_i. Note that a decrease of HR_i to lower values by age is also in concert with an increasing sympathetic dominance. Recent studies of Singh et al., using other nonlinear HRV measures, such as Approximate Entropy (ApEn) and Recurrence Quantification Analysis (RQA) analysed two age cohorts (”young” and ”elderly” subjects), also came to a similar conclusion [49], suggesting that ApEn and RQA could more sensitively measure age-dependent effects than linear parameters can do.

It is also worth emphasizing that similarly to the decrease of H_tp and H_i values, the centre of the VHF “cloud” also decreases as age progresses. This is a nontrivial correlation in the sense that, unlike the weight of the spectral components, the nonlinear HRV parameters do not directly measure the specific peaks considered a manifestation of ANS-related effects, but rather reflect the background noise structure, probably governed by coupled oscillations and fluctuations of other regulatory processes. Hence, they must carry some independent information, as well, to those HRV time- or frequency-domain parameters that measure the amplitude (power) of fluctuations.

4. Conclusions

Based on the intrinsic HR- and age-dependence of typical time- and frequency-domain HRV parameters and nonlinear HRV indices, we carried out a comparative study to reveal their fundamental connection, which is substantiated by basic mathematical considerations, and was experimentally demonstrated via the analysis of 24-hours ECG recordings of more than 200 healthy individuals.

We confirmed that the HR-dependence of one of the most frequently used time-domain parameters, RMSSD can be derived from Poincaré plots, as well, and that this “Master-Curve” representation seems to be superior to the conventional method. MC reflects both PNS and SNS effects, and represents the medium-term sympatho-vagal balance (from the hour- to the monthly scale). Momentary deviations from it (on the second to the minute scale) are indicative of the actual changes in, e.g., mood or stress [29], while long-term changes (on the yearly to the decade-scale) can be attributed to irreversible decline of the complex regulatory mechanism governing the dynamics of cardiac function (see also [27]).

In addition to this, we showed that the RMSSD(HR) and SDNN(HR) plots can also be reproduced from the Fourier-spectral components of the original RR time series, which is a natural manifestation of the Parseval theorem. As an illustrative tool for the visualization of the HRV spectra, we used a novel Fourier-map representation. By comparing the HR-dependence of the Fourier-map and the Master Curve, we found a new cloud-like structure of the HRV frequency spectrum. It is bunched mainly in the VHF band, and must be an important indicator of the effects of sympathetic nervous system on the heart rhythm. By an age-cohort evaluation, it was found that this “VHF band” is gradually shifted to lower HR regimes, as age progresses, and as such, it may be an important indicator of the age-related remodeling of the heart. We expect it to be utilized in the diagnostics of various cardiovascular diseases in the near future [45].

An anti-correlation between the DFA exponent (α) and sample entropy is apparent both on the individual and the cohort levels, which could be interpreted by frequency-spectral features of the HRV, as well. One can establish that these nonlinear indices report on the stochastic structure of HR fluctuations, contrary to the other time- or frequency-domain HRV measures that reflect direct PNS and SNS effects on the power of the HRV signal, hence, they carry independent information about the complex regulatory mechanism of the cardiac cycle. Nevertheless, the age-dependence of the extremum-places of these nonlinear measures seem to show correlation with the frequency-shift of the diffuse structure of the HRV spectrum, which can probably be assigned to increasing SNS dominance, as age progresses.

On the whole, we believe that our new findings regarding the phenomenology of HRV will shed more light on the context of various HRV measures, and will contribute to a deeper interpretation of the effects of external and internal factors associated to normal physiological or pathological phenomena.