The effect of age on the heart rate variability of healthy subjects

Abstract

In this work we study the characteristics of heart rate variability (HRV) as a function of age and gender. Our analysis covers a wider age range than that studied so far. It includes results previously reported in the literature and reveals behaviours not reported before. We can establish basic scale relationships in different HRV measurements. The mean value of the RR intervals shows a power-law behaviour independent of gender. Magnitudes such as the standard deviation or pNN50 show abrupt changes at around the age of 12 years, and above that age they show gender dependence, which mainly affects short-time (or high frequency) scales. We present a unified analysis for the calculation of the non-linear α and β parameters. Both parameters depend on age; they increase in the extremes of life and reach a minimum at around one year of age. These gender-independent changes occur at low frequencies and in scale ranges that depend on age. The results obtained in this work are discussed in terms of the effects of basal metabolic rate, hormonal regulation, and neuronal activity on heart rate variability. This work finally discusses how these findings influence the interpretation of HRV measurements from records of different lengths.

About the Authors:

Leopoldo Garavaglia

Roles Formal analysis, Writing – review & editing

Affiliation: Centro Científico Tecnológico (CCT) La Plata- CONICET, Instituto de Investigaciones Fisicoquímicas Teóricas y Aplicadas (INIFTA), Facultad de Ciencias Exactas, Universidad Nacional de La Plata, La Plata, República Argentina

Damián Gulich

Roles Formal analysis, Investigation, Methodology, Software, Writing – original draft, Writing – review & editing

Affiliation: Centro de Investigaciones Opticas, La Plata, República Argentina

Magdalena M. Defeo

Roles Data curation, Formal analysis, Investigation, Methodology, Resources, Validation, Writing – original draft, Writing – review & editing

Affiliation: Hospital Interzonal General de Agudos “Prof. Dr. Rodolfo Rossi”, La Plata, República Argentina

Julieta Thomas Mailland

Roles Data curation, Investigation, Methodology

Affiliation: Centro Científico Tecnológico (CCT) La Plata- CONICET, Instituto de Investigaciones Fisicoquímicas Teóricas y Aplicadas (INIFTA), Facultad de Ciencias Exactas, Universidad Nacional de La Plata, La Plata, República Argentina

Isabel M. Irurzun

Roles Conceptualization, Formal analysis, Funding acquisition, Project administration, Resources, Software, Supervision, Writing – original draft, Writing – review & editing

* E-mail: [email protected]

Affiliation: Centro Científico Tecnológico (CCT) La Plata- CONICET, Instituto de Investigaciones Fisicoquímicas Teóricas y Aplicadas (INIFTA), Facultad de Ciencias Exactas, Universidad Nacional de La Plata, La Plata, República Argentina

ORCID logo https://orcid.org/0000-0001-5357-2107

Introduction

Heart rate variability (HRV) is the physiological variation in the duration of cardiac cycles [1, 2]. Research into HRV began with the emergence of modern signal processing in the 1960s and 1970s, and has rapidly expanded in more recent times [3]. HRV is mainly controlled by the autonomic nervous system (ANS) through the interplay of sympathetic and parasympathetic neural activity especially at the sinus node [4, 5]. The autonomic nervous system interacts with receptors on the sinoatrial node (SAN) and the relative contributions of the ANS and SAN to HRV have been elucidated in [6, 7]. In general, HRV is influenced by many several factors such as chemical, hormonal and neural modulations, circadian changes, exercise, emotions, posture, and preload. The heart rate adaptation to changing factors is carried out by the activity of a variety of regulatory subsystems, and results in a complex linear and non-linear time behaviour, which changes with age and pathological conditions.

Several studies demonstrated age-related and gender-related variations in long-term HRV characteristics. It was reported that autonomic activities diminish with age in both genders and that gender-related variation in parasympathetic regulation decreases after the age of 50 years [8–14]. The hormonal influence on autonomic activity was analysed in [15–20]. Some investigators found that cardiac vagal activity declined during the luteal phase, and others reported that no changes in cardiovascular autonomic control occurred during the menstrual cycle of healthy fertile women. These studies were conducted on small populations and using HRV series too short to compare the results with ours.

HRV characteristics were proposed as predictors of the risk of premature mortality after myocardial infarction or development of congestive heart failure, also as diagnosis tools of autonomic dysfunction in diabetes, and as non-invasive tools for estimating the autonomic modulation of the cardiovascular system during stress, relaxation, or assessment of the effects of physical training on fitness level. For these reasons, the interest in HRV is growing both in clinical and physiological studies [21–28]. The importance of HRV as an index of the functional status of physiological control systems is currently recognized.

To compute HRV characteristics, many mathematical methods have been developed, which may be grouped into statistical, spectral, graphical, non-linear, complexity or information based [29–32].

In summary, the scientific and medical community has made a huge effort to obtain reliable measurements of the HRV characteristics in normal and pathological conditions. The study of the dependence of HRV on age and gender in a healthy population reveals the influence of factors such as body mass, basal metabolic rate, hormonal regulation, neuronal activity, etc. It also allows characterizing a control group to reference deviations associated with various pathologies.

There are many studies in the literature that together explore an age range of 1-99 years [2, 33–42].

In this paper we present a study that covers an age range of 0.08 (1 month)-99 years, increasing the age range by slightly more than a power of ten. The difference is significant for revealing scaling behaviours in the HRV measurements. In our analysis we include data obtained from [2, 33–35] that are representative of the age range 1-99 years and sufficient for the purpose of this work.

Finley et al. presented calculations of low frequency (LF), high frequency (HF), and total power indices for the range of 0.08-24 years [43]. LF and HF indices are not analysed in our work, and total power is equivalent to the standard deviation of HRV, whose dependence on age is discussed in this work. In that study, however, the data are grouped into age ranges, that is, age is not considered a continuous variable. Van den Berg et al. [44] presented a very large study including some of the indices that we calculate, but they examined 10-second electrocardiograms, which are too short to calculate scale factors by power and multifractal analyses, as we discuss in the last section.

We calculate indices in the time domain and scale factors (α and β) derived from both the multifractal and the power analyses, and a unified analysis of them is presented. There is a mathematical relationship between α and β, and our unified analysis suggests that the scale ranges depend on age. We also discuss α and β dependences on age in terms of the influence of neuronal activity on heart rate variability (due to either maturation of the autonomic nervous system or alterations in its interaction with the cells of the sinoatrial node).

Our work is limited to the analysis of the above-mentioned HRV characteristics [45, 46], but the methodological aspects reported in this work are general and should be considered to analyse any other measurements. By considering the HRV as a signal composed of multiple frequencies influenced by control mechanisms operating in different frequency ranges, there are two immediate aspects to consider: (i) the length of the HRV series is a variable coupled to the physiological variables, (ii) the HRV measurements must be determined in the same scale ranges in both the time and frequency domains. In other words, the length of the HRV series limits the scale range that can be analysed [47, 48]. There is a growing interest in using short and ultrashort HRV records because of their obvious clinical utility. To do this, it is necessary to understand the differences with 24-hour measurements, which have not yet been fully standardized. It is evident that by shortening the series, the frequencies under study are limited. A better understanding of the frequency composition of HRV and its dependence on age and gender can facilitate the use of short records.

This work is organised as follows: in Procedure we explain the methodology used, which is the same as that in [47, 48] and comparable to the one in [2, 33, 34]. The next section deals with the results; results of indices in the time domain are presented first, and then are results of scale factors α and β. The last section summarizes our conclusions.

Procedure

We analysed Holter recordings from healthy subjects who were recruited as volunteers after an exhaustive interview and clinical examination. Approval for this study was granted by the Ethics Committee of the National University of La Plata (UNLP) for data protection and privacy. In accordance with the Declaration of Helsinki, prior to participation in the study, all subjects were informed about the study and gave their written consent (in the case of children under 16, their parent’s or legal guardian’s consent was required).

The study included participants without clinical symptoms of disease, who were not on medication and whose electrocardiograms (ECG) were normal according to the criteria summarized in Table 1 [47, 48].

[Figure omitted. See PDF.]

Table 1. Criteria of normality for all Holter records obtained in the present work.

https://doi.org/10.1371/journal.pone.0255894.t001

Following these criteria, a total of 154 individuals were selected, who were aged between 1 month and 55 years (the age distribution is shown in Table 2).

[Figure omitted. See PDF.]

Table 2. Age distribution of the 154 individuals included in the research.

https://doi.org/10.1371/journal.pone.0255894.t002

Holter monitoring was recorded for 24 h with digital three-channel DMS300 7 and DMS300 3A recorders, and Galix recorders, using 3M electrodes [49]. The Galix recorders had a programmable read in sampling rate of 512 and 1024 Hz, and a write out sampling rate of 128 Hz. The DMS recorders had a sampling rate of 1024 Hz per channel for signal-averaged electrocardiography (SAEG) analysis, a read in sampling rate of 512 Hz, and a write out sampling rate of 128 Hz in the other cases. The signals were analysed with Galix software, and CardioScan 10.0, 11.0 software for DMS recorders. The error in the RR interval determination was of about 8 ms (twice the error in the determination of the R peak).

The records of cardiac events were automatically detected and classified by the Holter software, and then examined and corrected by two cardiologists. We employed the quality criteria established in [47, 48] for all time series used in the present work. Also stationarity was evaluated, and surrogate analysis was performed as in [50–53].

Time series of a total of 195 healthy participants were analysed. Of these time series, 154 corresponded to the participants mentioned above, 13 to [54, 55], and 28 to [54, 56]. Participants from [54, 56] were aged from 20 to 74 years (the age distribution is shown in Table 3), and 50% of them were females. Holter monitoring was recorded for 24 h with digital three-channel DMS recorders, the signals were analysed with CardioScan (so the error in the RR interval determination was also 8 ms), and then examined and corrected by two cardiologists. The HRV time series were selected after the evaluation according to quality criteria established in [47, 48]. For comparison we also analysed data from [2, 33–35]. The characteristics of the population studied in [2, 34] are shown in Tables 4 and 5 respectively. In [33] two hundred sixty healthy subjects 10 to 99 years old (112 males, 148 females) were recruited (the age distribution is shown in Table 6). Healthy subjects were defined as those without clinical evidence of organic disease in terms of medical history, physical examination, rest 12-lead electrocardiogram, routine blood chemistry profiles, and complete blood count. Oral contraceptives and non-steroidal anti-inflammatory agents were the only medications allowed. Subjects who exhibited abnormalities in one or more of the aforementioned categories were excluded from the study. Holter monitoring was recorded for 24h and analysed with different three recorders and software, and randomly selected Holter tapes were cross-analysed with the use of all programs. HRV and HR determinations from identical tapes were within 10% of each other. Recordings <20 h in duration, and those in which nonsinus beats comprised 10% of the total number of beats, were excluded.

[Figure omitted. See PDF.]

Table 3. Age distribution of the 41 individuals from [54–56].

https://doi.org/10.1371/journal.pone.0255894.t003

[Figure omitted. See PDF.]

Table 4. Clinical characteristics of the population used in [2].

https://doi.org/10.1371/journal.pone.0255894.t004

[Figure omitted. See PDF.]

Table 5. Clinical, lifestyle and laboratory characteristics of the population used in [34].

https://doi.org/10.1371/journal.pone.0255894.t005

[Figure omitted. See PDF.]

Table 6. Age and gender distribution of the individuals from [33].

https://doi.org/10.1371/journal.pone.0255894.t006

In total, we evaluated data from about 700 healthy subjects aged between 1 month and 99 years.

The age distributions in Tables 2, 3 and 6 are given as histograms dividing the population into subgroups. However, in our analysis we consider the age as a continuous variable and therefore we did not carry out an analysis of the size effect in these groups. The effect of size is considered in linear regressions as explained below.

The following linear indices in the time domain were calculated: the mean value of the RR intervals (<RR>), the standard deviation of the RR intervals (SDRR), the square root of the mean squared sum of differences of successive RR intervals (rMSSDRR), and the percentage of the intervals that vary by more than 50 ms from the previous interval (pNN50). These indices were calculated using the complete series, without detrending procedure. As we will discuss below, the length of the series affects its frequency composition (also does the detrending procedure) and this fact must be considered when comparing results obtained from short and long series. Our series are long enough so that the mentioned indices have values independent of the length of the series [47, 48].

As we mentioned in the Introduction, there are other widely used HRV indices, for example total power, and the LF and HF indices [33, 34, 43, 44]. The total power is equivalent to SDRR. The LF and HF indices are calculated over specific frequency ranges. One of our conclusions points to the need to improve the definition of the frequency ranges studied, which would be dependent on age. For this reason we focus on the scale factors and the relationship between them, as explained below.

It is known that HRV follows a power-law behaviour in the frequency domain (1/beat), which is manifested in the power spectrum dependence on frequency as(1)

The values for β were determined in this work by averaging the power spectra of successive time series segments, each consisting of 4096 datapoints. This procedure, which allows eliminating high frequency fluctuations, is detailed in [47]. Power spectra were calculated using fast Fourier transform (FFT) after normalizing the series to zero mean value. The method to calculate β used the optimal range for each participant (roughly in the range 2 × 10−4–5 × 10−2Hz) [47]. No detrending procedure was used. In the next section we compare our results with those obtained in [34]. In [34]the slope (β) was calculated from the frequency range of 10−4 to 10−2 Hz from the 24-hour ECG recordings. The point power spectrum was logarithmically smoothed in the frequency domain, and the power was integrated into bins spaced 0.0167 log(Hz) apart. A line-fitting algorithm of log(power) on log(frequency) was then applied to the power spectrum between 10−4 and 10−2 Hz, and the slope of this line was calculated.

The detrended fluctuation analysis (DFA) is widely used to characterize the fractal dynamics of a system from which a time series has been measured [57, 58], and it is the most popular approach to detect the presence of long-term memory in data [59]. By the DFA technique a scale exponent α is determined, which is related to β as: β = 2α−1. However, we rarely find a verification of this relationship in scientific papers. There are several reasons for this omission. The value of α depends on the scale ranges chosen to calculate it, and there is no agreement on such ranges. The scale ranges for calculating β should be compatible with those used for α, but the reliable calculation of the power spectrum usually requires more data.

The DFA follows five steps [60–62]

* Step 1: given a time series S = {xt, t = 1, …, N}, with N being the number of equidistant observations, the cumulative sum of the data , with k = 1, …, N and , is considered.

* Step 2: this profile is divided into ⎣N/s⎦ non-overlapping windows of equal length s (⎣a⎦ denotes the largest integer less than or equal to a).

* Step 3: a local polynomial fit yν,m(k) of degree m is fitted to the profile for each window ν = 1, …, ⎣N/s⎦. The degree of the polynomial can be varied to eliminate constant (m = 0), linear (m = 1), quadratic (m = 2) or higher order trends of the profile. It is customary to indicate the degree of detrending by including it in the title of the technique (DFA- m).

* Step 4: the variance of the detrended time series is evaluated by averaging over all datapoints k in each segment ν, , for ν = 1, …, ⎣N/s⎦.

* Step 5: the DFA fluctuation function is obtained by averaging over all segments and taking the square root, .

By repeating this procedure for different values of s, the s-dependence of Fm s is obtained.

If the time series has long-range power-law correlations, Fm(s) scales as(2)for a certain range of s [62]. The scaling exponent α is estimated by the slope of the best linear regression in a double logarithmic plot. The long-range correlations embedded in the time series are quantified by this exponent: if α > 1/2, consecutive increments tend to have the same sign, so the processes are persistent. If α < 1/2, consecutive increments are more likely to have opposite signs, and it is said that the processes are anti-persistent. For uncorrelated data α = 1/2 is obtained [63].

We calculated the fluctuation functions F2(s) with detrending degrees m = 1 − 4 (scales s ranging from 5 to 200 datapoints). For m = 3, 4, detrending artifacts appear, and therefore m values should be limited to 1 or 2 for a correct estimation of α. To compare with previous results, we used three sets of scales s to fit all fluctuation functions, thus finding three scaling exponents (Eq (2)):

* α0 for scales 5 ≤ s ≤ 10 (very high frequency exponent) [34]

* α1 for scales 10 ≤ s ≤ 50 (high frequency exponent) [35]

* α2 for scales 50 ≤ s ≤ 200 (medium frequency exponent) [35]

Note that the scale ranges used to calculate αi(i = 0, 1, 2) are somewhat different from those used in Peng’s original article [58]. However, these differences do not imply significant changes in the values of αi or modify the conclusions of this work. We used the scales that were employed in the papers whose data were used for comparison. In fact, in this work we show that the scale ranges are dependent on age and must be determined consistently in both the time and frequency domains.

In Results the relationships of the linear indices with age are presented. Age is considered as a continuous variable. A logarithmic transformation is then applied to both the independent and dependent variables to reveal scale relationships between the variables, which is manifested as a straight line in a log-log plot. The scale factors are determined by performing a linear regression fit using the least squares method, and the significance of the results is evaluated by an analysis of variance (ANOVA). The size effects can be evaluated either through Cohen’s d index or the regression coefficient (R). We used R in this work, and the values obtained indicate that the size effects are negligible in the relationships proposed in this work (Eqs 3–6). Other adjustments previously proposed in the scientific literature were also evaluated, and the ANOVA tables showed that the scale relationships were the most acceptable. The revelation of the scale relationships in this work is a consequence of the widening of the age range and the consideration of age as a continuous variable. When analysing the influence of gender, we performed a covariance analysis (ANCOVA). The size effects were also negligible in these cases.

The analyses of the linear indices included data from other studies that show differences in the methodology or in the population studied. As far as we could verify, these differences did not significantly affect our results (see Results).

Finally, although age is considered a continuous variable, our study is cross-sectional. In the proposed adjustments for the linear indices we then included the prediction intervals (see Fig 1). A prediction interval is an estimate of a range of values in which a future observation will occur with a certain probability, given what has already been observed. Thus, in a hypothetical longitudinal study, the indices of a healthy individual should evolve within the prediction intervals.

[Figure omitted. See PDF.]

Fig 1. Age dependence of some time-domain indices: (A)<RR>, (B) SDRR, (C) rMSSDRR, (D) pNN50 Filled circles correspond to data obtained in this work, open squares to data from [2], open circles to data from [33], and open triangles to data from [34].

Dashed lines indicate the 95% prediction intervals.

https://doi.org/10.1371/journal.pone.0255894.g001

The dependence of the α and β scale factors on age is presented in linear-log plots. In these cases, growth trends and the existence of critical points (maximum or minimum) are analysed, and the application of a logarithmic transformation to age was done in order to better visualize the behaviour at small age values. We do not propose a specific functional form for α and β dependence on age in this work. The prediction intervals could be roughly calculated by locally approximating the function by line segments. We did this for β in Fig 2 and α2 in Fig 3 for the reasons we discuss in Results.

[Figure omitted. See PDF.]

Fig 2. β dependence on age.

Filled circles correspond to data obtained in this work, open circles to data from [34]. Dashed lines indicate the approximate 95% prediction interval (see Procedure).

https://doi.org/10.1371/journal.pone.0255894.g002

[Figure omitted. See PDF.]

Fig 3. α0 (A), α1(B), and α2(C) dependence on age.

Filled circles correspond to data obtained in this work, and open squares and open circles correspond to data from [34], for α0 and α1 respectively. Note that in [34]α1 was determined for s > 11. In (C) the dashed lines indicate the approximate 95% prediction interval (see Procedure).

https://doi.org/10.1371/journal.pone.0255894.g003

Results

Time-domain indices

Fig 1 shows the dependence of time domain indices on age, including data from [2, 33, 34].

Data show good agreement among these indices, validating the general treatment of the measurements. Our data cover a wider age range, revealing novel tendencies. While <RR> exhibits a monotonic behaviour, SDRR, rMSSDRR, and pNN50 show an abrupt change at about the age of 12 years. In Fig 1 log-log scales reveal the existence of scale laws, which appear as linear behaviours after a logarithmic transformation. Results have been fitted to scaling laws as follows:(3) (4) (5) (6)where x is the age in years, and the statistics of the adjustments are shown in Table 7.

[Figure omitted. See PDF.]

Table 7. Statistics of Eqs 3, 4, 5 and 6.

https://doi.org/10.1371/journal.pone.0255894.t007

The threshold x = 12 was determined for SDRR, and the same value was assumed valid for rMSSDRRand pNN50. We made linear fits by changing the cut-off value, and x = 12 was chosen because the linear fits on both sides were the best. Though the power-law adjustments for ages above 12 years are statistically worse than those for ages below 12 years, they are still better than or equal to other linear or quadratic fits performed on the same data sets.

Quadratic fits over the complete data set have been tested, but they were statistically less reliable than those presented in this work. Table 8 shows the analysis of variance (ANOVA) tables of different fits performed in this work.

[Figure omitted. See PDF.]

Table 8. ANOVA tables for different fits over the complete data set tested in this work.

https://doi.org/10.1371/journal.pone.0255894.t008

We find slight differences by gender only for ages over 12. Table 9 summarizes the power-law parameters for each case. We performed an analysis of covariance (ANCOVA) to assess the influence of gender and the p-values of the covariate are shown in the last column.

[Figure omitted. See PDF.]

Table 9. Power law adjustments by gender for participants over 12 years old.

https://doi.org/10.1371/journal.pone.0255894.t009

Note that the R values do not improve when discriminating the population by gender. We also made adjustments for x > 12 only for the population in [33, 35], and the values of R were similar. A similar result was reported in [34]. It is evident that there is great individual variability in HRV indices in the population over 12, which we discuss in the last section.

Scale factors

Fig 2 shows the β dependence on age in a linear-log plot in order to better visualize the behaviour at small age values. For comparison we included data from [34]. The agreement between the two data sets is reasonable, but our data cover a wider age range than that reported so far. A novel non-monotonic behaviour is observed with an increase of β values in the extremes of life.

Fig 3 shows the dependence of α0, α1, and α2 (m = 1) on age, including data from [34].

Again the agreement between the two data sets is reasonable. In addition, the exponents exhibit a non-monotonous behaviour. The values of α0 are widely dispersed, and their calculation is discussed below;α1 and α2 have a minimum value at around one year of age.

Data in Ref [35] can be compared with α0 and α1. Though the general behaviour with age is similar, these data show more dispersion and decrease more rapidly with age, reaching values close to 0.5 for young adults. This type of deviation was observed in our data with the increase of m. High frequency scale factors are very sensitive to the increase of m.

We also analysed the relationship between α and β. Of the three scale factors, the one that correlates best with beta is α2. In Fig 4 we plot the quotient q as a function of age, being(7)

[Figure omitted. See PDF.]

Fig 4. dependence on age.

https://doi.org/10.1371/journal.pone.0255894.g004

If the scale relationship between α2 and β is fulfilled, q should be equal to 1, but the values of q deviate as age decreases. We discuss the possible causes of this discrepancy in the next section.

Discussion and conclusions

In this work we analyse data from 24-hour Holter records for a total of 700 healthy subjects, 154 of whom were volunteers recruited for this study. The analysis covers the largest age range analysed so far. As we mentioned before, many other works in the literature consider HRV dependence on age and gender, covering an age range of 1-99 years [2, 33–44]. Our work increases the age range by slightly more than a power of ten—a significant difference for detecting scaling behaviour. As we mentioned in the Introduction, Finley et al. [43] also reported data for ages as early as 1 month. Of the indices calculated in their work, total power is comparable to HRV standard deviation (see ii). It should be noted that in [43] the records analysed have a length of 4 h during sleep and 1 h during wakefulness. According to the quality criteria established in [47, 48], those series are too short to ensure the robustness of the calculations (that is, to produce results independent of the length of the series). Only the longest period could then be used for a qualitative comparison with our work. Another important methodological difference is the report of the results grouped by age ranges, that is, the consideration of age as a discrete variable instead of a continuous one. Van den Berg et al. [44] also presented a very large study that included an age range from 1 month, but it was performed on 10-second electrocardiograms, which are too short to allow calculating α and β indices. The use of short records is undoubtedly desirable for their practical utility, but the length of the series also affects their frequency composition, and therefore the comparison must be careful, as discussed below.

We calculate indices in the time domain: <RR>, SDRR, rMSSDRR, and pNN50, measurements in the frequency domain (β) and multifractal scale fractals (α).

Our work is limited to the analysis of the above-mentioned HRV measurements; however, the methodological aspects presented in this work are useful to analyse any other measurements.

Our main conclusions are the following:

1. (i)<RR> follows a scaling relationship with age, which is independent of gender. This power-law behaviour reveals a self-similar growth process, which can be related to the basal metabolic rate (BMR) [64–70]. Allometric relationships for <RR> with body mass (BM) or the BMR through Kleiber’s law can be written as:(8) (9)where it has been determined that a = 1/4 for mammalians, and the value of ddepends on the scale factor in the Kleiber’s law [65]. The d value was initially established equal to 3/4, though it has recently been questioned based on empirical evidence that showed a wealth of robust exponents deviating from 3/4 [69, 70]. In Eq (3) the scale relationship is established with age with a scale factor of ∼1/8.

An additional important aspect is that this scaling relationship remains unchanged throughout life. One should expect deviations at an even earlier age, and it would therefore be interesting to analyse the validity of this relationship from conception and not from birth.

1. (ii) Statistical measures such as SDRR, rMSSDRR, and pNN50 show an abrupt change at the age of 12 years. For simplicity we assume the same cut-off (independent of gender) for all measures, but this assumption should be further studied. The threshold x = 12 was then determined for SDRR looking for the best the linear fits on both sides of the threshold. Note in Fig 1 that data reported in [34] already indicated the existence of a maximum. But to attribute this maximum to a change in a scale behaviour requires analysing wide age ranges and using age as a continuous variable, as was done in this work. Quadratic fits over the complete data have been tested, but they were statistically less reliable than those presented in this work. Regarding data reported for HRV during sleep in [43], Finley et al. distinguish between quiet and active states. The maximum is clearly observed in quiet sleep between the ages of 6 and 11, while in the active state there is a continuous decrease. To compare with our work, it is appropriate to consider the two states together, and the maximum then occurs near the age of 11 years. However, a closer comparison is not possible due to the methodological differences mentioned above.

Changes of SDRR, rMSSDRR, and pNN50 at the age of 12 years are related to HRV changes at high frequencies, i.e. to alterations in the frequency composition of HRV. We note that the scale law of <RR> does not change, while those of the other measures do, i.e. over time there is a redistribution at lower frequencies, without altering the mean value of the series.

An increase in individual variation is also observed above 12 years of age, which continues even when discriminating by gender. As we mentioned in the results, this large individual variation, which is manifested in low R values, cannot be attributed to the differences between the studies used in this work because they were also observed in each of these studies. On the other hand, a large interindividual variation of autonomic function has also been described in healthy subjects [71–74]. If a process of redistribution does indeed occur in the frequency composition of HRV with age, then individual changes in the onset and progression of this process may explain the observed variability. In [6, 7] the authors studied the influence of both the ANS and the SAN receptors on HRV. They found that the deterioration of mechanisms intrinsic to the SAN leads to an increased <RR> and HRV. The extrinsic mechanisms (ANS input) compensate for deteriorated intrinsic mechanisms to preserve <RR>, but this compensation is associated with a reduction in the HRV.

Below the age of 12, results are gender-independent, while above that age, there seems to be a slight dependence on gender. The results suggest that the frequency composition of HRV is also altered. In fact, this alteration seems to occur at high frequencies because β values, which are determined at low frequencies, do not depend on gender (see below). Studies in [16] reported sex- differences in the HRV under stress. The authors suggested that men and women have different autonomic “strategies” for dealing with stress, such that men rely more on the hypothalamic-pituitary-adrenal (HPA) axis activation, and women on parasympathetic withdrawal. The researchers did not report differences in the baseline state (without stress), but these differences could not be detected due to the age range (18-60 years) and the length of the records used. We believe that we are measuring a baseline state, the one that was not detected in [16].

There are studies that support a potential protective role of estradiol in the parasympathetic cardiac tone [17, 18]. The effect of progesterone and its potent metabolite allopregnanolone on ANS is uncertain. Our study was not controlled for menstrual cycle phase. We think that the differences found in this work are beyond those that can be attributed to the cycle by the age range and the magnitude of the population analysed, and by the diversity of the origins of the information. In other words, the effects of the cycle phases should be averaged.

Finally, if the length of the electrocardiographic records limits the frequency range that can be analysed, and this range is also influenced by gender, it is clear that long and short records will differ in the information they provide and also in their clinical utility. Even though short records may not necessarily replace long ones in all applications, additional understanding of the information provided by each type would be necessary.

1. (iii) We performed a unified analysis of the scale factors (α and β) derived from both the multifractal and the spectral power analyses. Because of the mathematical relationship between these factors, the calculation was done by two independent methods. Whereas for the calculation of α we respected the scale ranges used in the literature, the method developed to calculate β used the optimal range for each participant. Of the three scale factors derived from the multifractal analysis, the one that correlates best with β is α2. Although the relationship between α and β can be then verified reasonably well, at least for people older than 5, there is a systematic deviation at an earlier age because the scale ranges are age-dependent. But more importantly, as a methodological change, the scaling hypothesis must be tested in the same range in both the time and the frequency domains. Therefore the determination of the scale factors must be based on simultaneous and independent calculations of α and β. The scale factors α and β exhibit a non-monotonic relationship with age, increasing in the extremes of life and reaching a minimum at around one year of age. This behaviour, evident in α1 and α2, has not been reported before. The values of α0 are widely dispersed due to the narrow scale range (5 < s < 10) used in the calculation. In this scale range, a certain dependence on gender should be evidenced, but the dispersion of the results prevents noticing it. We attribute the changes in α2 and β to the neuronal control on the HRV. As they are calculated at low frequencies (long-range scale), hormonal effects are negligible, and results are gender-independent. In [47, 48] we showed that an increase in β values can be attributed to a predominance of the sympathetic system over the parasympathetic one. At early ages, this predominance is due to the lack of maturation of the parasympathetic modulation, while at older ages it is explained by a decrease in the parasympathetic modulations.

Our knowledge of the population characteristics is limited, we do not have information on the menstrual status of all women, the socioeconomic status of the population, the BMI, etc. All studies met minimum requirements, patients were examined by clinicians and minimum quality requirements were established on Holter records and time series. However, the clinical normality criteria were left to the intervening physicians and our study was based solely on the analysis of Holter records. Our results do not allow us to distinguish the contributions of ANS and SAN [6, 7]. In this context, the sympathetic and parasympathetic modulations mentioned above are those detected at the output of the SAN.

In summary, we can establish basic scale relationships in different HRV measurements performed on long records. Above the age of 12, there is a gender dependence that mainly alters the high frequency region. Changes attributable to neuronal activity (including SAN) are observed at low frequencies. These differences must be considered when interpreting measurements based on recordings of a few seconds.

The scale factors should be calculated individually in age-dependent scale ranges and equally in both time and frequency domains.

In [75] we presented the basis for the development of implantable devices that mimic RR variability of the healthy heart. The present work has direct implications for that project that is now in progress, because that development requires a comprehensive characterization of HRV, such as the one this study presents.

Citation: Garavaglia L, Gulich D, Defeo MM, Thomas Mailland J, Irurzun IM (2021) The effect of age on the heart rate variability of healthy subjects. PLoS ONE 16(10): e0255894. https://doi.org/10.1371/journal.pone.0255894

References

1. Voss A, Schroeder R, Heitmann A, Peters A, Perz S. Short-Term Heart Rate Variability-Influence of Gender and Age in Healthy Subjects. PLoS ONE 2015; 10(3): e0118308. pmid:25822720

2. Bobkowski W, Stefaniak ME, Krauze T, Gendera K, Wykretowicz A, Piskorski J, et al. Measures of heart rate variability in 24-h ECGs depend on age but not gender of healthy children. Front. Physiol 2017; 8:311. pmid:28572771

3. McCraty R, Shaffer F Heart rate variability: new perspectives on physiological mechanisms, assessment of self-regulatory capacity, and health risk. Global Advances in Health and Medicine 2015; 4(1):46–61. pmid:25694852

4. Task Force. Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology. Circulation. 1996; 93(5): 1043–1065.

5. Schwab JO, Eichner G, Schmitt H, Weber S, Coch M, Waldecker B The relative contribution of the sinus and AV node to heart rate variability. Heart. 2003; 89(3): 337–338. pmid:12591853

6. Rosenberg AA, Weiser-Bitoun I, Billman GE, Yaniv Y Signatures of the autonomic nervous system and the heart-’s pacemaker cells in canine electrocardiograms and their applications to humans. Scientific Reports 2020; 10:9971. pmid:32561798

7. Yaniv Y, Ahmet I, Tsutsui K, Behar J, Moen JM, Okamoto Y, et al. Deterioration of autonomic neuronal receptor signaling and mechanisms intrinsic to heart pacemaker cells contribute to age-associated alterations in heart rate variability in vivo. Aging Cell 2016; 15: 716–724. pmid:27168363

8. Jensen-Urstad K, Storck N, Bouvier F, Ericson M, Lindblad LE, Jensen-Urstad M Heart rate variability in healthy subjects is related to age and gender. Acta Physiol Scand. 1997; 160(3): 235–241. pmid:9246386

9. Pikkujamsa SM, Makikallio TH, Sourander LB, Raiha IJ, Puukka P, Skytta J, et al. Cardiac Interbeat Interval Dynamics From Childhood to Senescence. Circulation 1999; 100(4):393–399. pmid:10421600

10. Goldberger AL Non-linear dynamics for clinicians: chaos theory, fractals, and complexity at the bedside. Lancet. 1996; 347(9011): 1312–1314. pmid:8622511

11. Gribbin B, Pickering TG, Sleight P, Peto R Effect of age and high blood pressure on baroreflex sensitivity in man. Circ Res. 1971; 29(4): 424–431. pmid:5110922

12. Voss A, Schulz S, Schroeder R, Baumert M, Caminal P Methods derived from nonlinear dynamics for analysing heart rate variability. Philos Transact A Math Phys Eng Sci. 2009; 367: 277–296. pmid:18977726

13. Kleiger RE, Stein PK, Bigger JT Jr. Heart rate variability: measurement and clinical utility. Ann Noninvasive Electrocardiol. 2005; 10(1): 88–101. pmid:15649244

14. Porta A, Guzzetti S, Montano N, Furlan R, Pagani M, Malliani A, et al. Entropy, entropy rate, and pattern classification as tools to typify complexity in short heart period variability series. IEEE Trans Biomed Eng. 2001; 48(11): 1282–1291. pmid:11686627

15. Koifman R, Dayan L, Ablin JN, Jacob G Cardiovascular Autonomic Profile in Women With Premenstrual Syndrome, Frontiers in Physiology 2018; 9:1384. pmid:30327616

16. Hamidovica A, Van Hedgerb K, Choia SH, Flowersa S, Wardlec M, Childs E Quantitative meta-analysis of heart rate variability finds reduced parasympathetic cardiac tone in women compared to men during laboratory-based social stress, Neuroscience and Biobehavioral Reviews 2020; 114: 194–200.

17. Bai X, Li J, Zhou L, Li X Influence of the menstrual cycle on nonlinear properties of heart rate variability in young women. Am. J. Physiol.-Heart Circul. Physiol. 2009; 297: H765–H774. pmid:19465541

18. McKinley P, King A, Shapiro P, Slavov I, Fang Y, Chen I, Jamner L, et al. The impact of menstrual cycle phase on cardiac autonomic regulation. Psychophysiology 2009; 46: 904–911. pmid:19386049

19. Tenan MS, Brothers RM, Tweedell AJ, Hackney AC, Griffin L. Changes in resting heart rate variability across the menstrual cycle. Psychophysiology 2014; 51:996–1004. pmid:24942292

20. von Holzen JJ, Capaldo G, Wilhelm M, Stute P. Impact of endo- and exogenous estrogens on heart rate variability in women: a review. Climacteric 2016; 19: 222–228. pmid:26872538

21. Bigger JT, Fleiss JL, Steinman RC, Rolnitzky LM, Kleiger RE, Rottman JN. Frequency domain measures of heart period variability and mortality after myocardial infarction. Circulation 1992; 85(1): 164–171. pmid:1728446

22. La Rovere MT, Bigger JT Jr, Marcus FI, Mortara A, Schwartz PJ. Baroreflex sensitivity and heart rate variability in prediction of total cardiac mortality after myocardial infarction. ATRAMI (Autonomic Tone and Reflex After Myocardial Infarction) Investigators. Lancet 1998; 351(9101):478–484. pmid:9482439

23. Stein PK, Reddy A. Non-linear heart rate variability and risk stratification in cardiovascular disease. Ind. Pacing Electrophysiol. J. 2005; 5(3): 210–220. pmid:16943869

24. Guzik P, Piskorski J, Barthel P, Bauer A, Muller A, Junk N, et al. Heart rate deceleration runs for postinfarction risk prediction. J. Electrocardiol. 2012; 45(1): 70–76. pmid:21924431

25. Guzik P, Piskorski J, Krauze T, Schneider R, Wesseling KH, Wykretowicz A, et al. Correlations between the Poincaré plot and conventional heart rate variability parameters assessed during paced breathing. J. Physiol. Sci. 2007; 57(1): 63–71. pmid:17266795

26. Guzik P, Piskorski J, Krauze T, Wykretowicz A, Wysocki H. Partitioning total heart rate variability. Int. J. Cardiol. 2010; 144(1): 138–139. pmid:19178958

27. Patel VN, Pierce BR, Bodapati RK, Brown DL, Ives DG, Stein PK. Association of Holter-derived heart rate variability parameters with the development of congestive heart failure in the cardiovascular health study. JACC Heart Fail. 2017; 5(6): 423–431. pmid:28396041

28. Akinci A, Celiker A, Baykal E, Tezic T. Heart rate variability in diabetic children: sensitivity of the time and frequency domain methods. Pediatr. Cardiol. 1993; 14(3): 140–146. pmid:8415215

29. Huikuri HV, Makikallio TH, Perkiomaki J. Measurement of heart rate variability by methods based on nonlinear dynamics. J. Electrocardiol. 2003; 36 Suppl: 95–99. pmid:14716599

30. Goldberger AL, Amaral LA, Hausdorff JM, Ivanov P, Peng CK, Stanley H.E. Fractal dynamics in physiology: alterations with disease and aging. Proc Natl Acad Sci U S A. 2002; 99 Suppl 1: 2466–2472. pmid:11875196

31. Huikuri HV, Perkiomaki JS, Maestri R, Pinna GD. Clinical impact of evaluation of cardiovascular control by novel methods of heart rate dynamics. Philosophical transactions Series A, Mathematical, physical, and engineering sciences. 2009; 367(1892): 1223–1238. pmid:19324705

32. Rajendra Acharya U, Joseph KP, Kannathal N, Lim CM, Suri JS. Heart rate variability: a review. Med Biol Eng Comput. 2006; 44(12): 1031–1051. pmid:17111118

33. Umetani K, Singer DH, McCraty R, Atkinson M. Twenty-four hour time domain heart rate variability and heart rate: relations to age and gender over nine decades. JACC 1998; 31(3): 593–601. pmid:9502641

34. Pikkujamsa S. Heart Rate Variability and Baroreflex Sensitivity in Subjects without Heart Disease. Effects of age, sex and cardiovascular risk factors. PhD Thesis. Oulu University. Finland. 1999. Available from http://jultika.oulu.fi/files/isbn9514252276.pdf

35. Shiogai Y, Stefanovska A, McClintock PVE. Nonlinear dynamics of cardiovascular ageing. Physics Reports 2010; 488 (2-3): 51110. pmid:20396667

36. Shaffer F, Ginsberg JP. An Overview of heart rate variability metrics and norms, Front. Public Health 2017; 5:258.

37. Silvetti MS, Drago F, Ragonese P. Heart rate variability in healthy children and adolescents is partially related to age and gender. International Journal of Cardiology. 2001; 81(2-3):169–174. pmid:11744133

38. Tegegne BS, Man T, van Roon AM, Riese H, Snieder H. Determinants of heart rate variability in the general population: the lifelines cohort study. Heart Rhythm. 2018; 15(10):1552–1558. pmid:29753022

39. Plaza -Florido A, Alcantara JMA, Migueles JH, Amaro -Gahete FJ, Acosta FM, Mora -Gonzalez J, et al. Inter and intra -researcher reproducibility of heart rate variability parameters in three human cohorts. Scientific Reports 2020; 10: 11399, https://doi.org/10.1038/s41598-020-68197-7 pmid:32647148

40. Bonnemeier H, Wiegand UKH., Brandes A, Kluge N, Katus HA, Richardt G, et al. Circadian profile of cardiac autonomic nervous modulation in healthy subjects: differing effects of aging and gender on heart rate variability. Journal of Cardiovascular Electrophysiology 2003; 14 (8): 791–799. pmid:12890036

41. Finley JP, Nugent ST, Hellenbrand W. Heart-rate variability in children. Spectral analysis of developmental changes between 5 and 24 years. Can. J. Physiol. Pharmacol 1987; 65(10): 2048–2052. pmid:3427543

42. Reardon M, Malik M. Changes in heart rate variability with age. Pacing and Clinical Electrophysiology, 1996; 19(11): 1863–1866. pmid:8945057

43. Finley JP, Nugent ST. Heart rate variability in infants, children and young adults. Journal of the Autonomic Nervous System 1995; 51(2):103–108. pmid:7738283

44. van den Berg ME, Rijnbeek PR, Niemeijer MN, Hofman A, van Herpen G, Bots ML, et al. Normal values of corrected heart rate variability in 10-second electrocardiograms for all ages. Front. Physiol. 2018; 9:424. pmid:29755366

45. Zaza A, Lombardi F. Autonomic indexes based on the analysis of heart rate variability: a view from the sinus node. Cardiovascular Research 2001, 50: 434–442. pmid:11376619

46. Costa MD, Goldberger AL. Heart rate fragmentation: using cardiac pacemaker dynamics to probe the pace of biological aging. Am J Physiol Heart Circ Physiol 2019; 316: H1341–H1344. pmid:30951362

47. Andrés DS, Irurzun IM, Mitelman J, Mola EE. Increase in the embedding dimension in the heart rate variability associated with left ventricular abnormalities. Appl. Phys. Lett. 2006; 89, 144111

48. Irurzun IM, Mola EE. Heart Rate Variability: A View from Chaos Theory. 1st. Edition. LAP Lambert Academic Publishing, 2015.

49. ACC/AHA Guidelines for Ambulatory Electrocardiography, Journal of the American College of Cardiology. 1999; 34(3): 912–948. pmid:10483977

50. Theiler J, Eubank S, Longtin A, Galdrikian B, Farmer JD. Testing for nonlinearity in time series: The method of surrogate data. Physica D 1992; 58(1-4): 77–94.

51. Schreiber T, Schmitz A. Improved surrogate data for nonlinearity tests. Phys. Rev. Lett. 1996; 77(4): 635–638. pmid:10062864

52. Kugiumtzis D. Test your surrogate data before you test for nonlinearity. Phys. Rev. E. 1999; 60 (3): 2808–2816. pmid:11970086

53. Popivanov D, Mineva A. Testing procedures for non-stationarity and nonlinearity in physiological signals. Mathematical Biosciences 1999; 157(1-2):303–320. pmid:10194935

54. Goldberger A, Amaral L, Glass L, Hausdorff J, Ivanov P.C, Mark R, et al. PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals. Circulation. 2000; 101 (23): e215–e220. pmid:10851218

55. Normal Sinus Rhythm RR Interval Database.

56. MIT-BIH Normal Sinus Rhythm Database.

57. Peng CK, Buldyrev SV, Havlin S, Simons M, Stanley HE, Goldberger AL. Mosaic organization of DNA nucleotides. Phys. Rev. E 1994; 49(2):1685–1689. pmid:9961383

58. Peng CK, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationaty heartbeat time series. Chaos 1995; 5(1): 82–87. pmid:11538314

59. Grech D, Mazur Z. On the scaling ranges of detrended fluctuation analysis for long-term memory correlated short series of data. Physica A 2013; 392(10): 2384–2397.

60. Kantelhardt JW, Koscielny-Bunde E, Rego HHA, Havlin S, Bunde A. Detecting long-range correlations with detrended fluctuation analysis. Physica A 2001; 295(3-4): 441–454.

61. Kantelhardt JW, Zschiegner SA, Koscielny-Bunde E, Havlin S, Bunde A, Stanley HE. Multifractal detrended fluctuation analysis of nonstationary time series. Physica A 2002; 316(1-4): 87–114.

62. Gulich D, Zunino L. A criterion for the determination of optimal scaling ranges in DFA and MF-DFA. Physica A 2014; 397: 17–30.

63. Feder J, Fractals, Plenun Press, New York, 1988.

64. Bassil G, Zarzoso M, Noujaim SF. Allometric scaling of electrical excitation and propagation in the mammalian heart. J Theor Biol. 2017; 419: 238–242. pmid:27686596

65. Kleiber M. Body size and metabolic rate. Physiol Rev. 1947; 27(4): 511–541. pmid:20267758

66. Kolokotrones T, Savage V, Deeds EJ, Fontana W. Curvature in metabolic scaling. Nature. 2010;464(7289):753–756. pmid:20360740

67. Noujaim SF, Lucca E, Muñoz V, Persaud D, Berenfeld O, Meijler FL, et al. From mouse to whale: a universal scaling relation for the PR Interval of the electrocardiogram of mammals. Circulation. 2004;110(18):2802–2808.

68. Noujaim SF, Berenfeld O, Kalifa J, Cerrone M, Nanthakumar K, Atienza F, et al. Universal scaling law of electrical turbulence in the mammalian heart. Proc Nat Acad Sci USA. 2007;104(52):20985–20989. pmid:18093948

69. Ballesteros FJ, Martinez VJ, Luque B, Lacasa L, Valor E, Moya A On the thermodynamic origin of metabolic scaling. Scientific Reports 2018; 8:1448. pmid:29362491

70. Sánchez S, Barragán J. Metabolically active weight: between Kleiber’s law and the second law of thermodinamics Rev. Argent. Endocrinol. Metab. 2011; 48(3): 136–142.

71. Abdel-Rahman AR, Merrill RH, Wooles WR. Gender-related differences in the baroreceptor reflex control of heart rate in normotensive humans. J Appl Physiol 1994; 77: 606–613. pmid:8002506

72. Kupari M, Virolainen J, Koskinen P, Tikkanen MJ. Short-term heart rate variability and factors modifying the risk of coronary artery disease in a population sample. Am J Cardiol 1993 72: 897–903. pmid:8213546

73. Huikuri HV, Kessler KM, Terracall E, Castellanos AX, Linnaluoto M, Myerburg RJ. Reproducibility and circadian rhythm of heart rate variability in healthy subjects. Am J Cardiol 1990 65: 391–393. pmid:2301268

74. Molgaard H, Hermansen K, Bjerregaard P. Spectral components of short-term RR interval variability in healthy subjects and effects of risk factors. Eur Heart J 1994 15: 1174–1183. pmid:7982416

75. Defeo MM. Influencia del modo de estimulación y de la tasa de captura sobre las propiedades complejas de la variabilidad de la frecuencia cardiaca en pacientes con marcapasos implantables. PhD Thesis. University of La Plata. Argentina. 2021.