Increase in random component of heart rate variability coinciding with developmental and degenerative stages of life

Processes of biological adaptation may be considered as attractor selections (Kashiwagi et al 2006). In order to adapt to a new environment or situation, organisms seem to explore the most convenient state while moving between attractors by random fluctuation (Yanagida et al 2007). For example, children find behaviors that are more suitable for situations through many random actions. If this is true for biological adaptions in general, physiological signals, particularly those reflecting brain activities, may increase the random components at the stages of life where the needs for adaption increase, i.e. during early developmental and final degenerative stages.

To examine this hypothesis, we measured the age-dependent changes in the random component in heart rate variability (HRV) and compared them with the values in patients with a developmental disorder (DD). Because HRV is known to disappear with complete cardiac autonomic denervation (Sands et al 1989), it is thought to originate from the brain and to reflect its activities to regulate the heart through the autonomic nervous system (Camm et al 1996, Hayano 2016). To assess the quantity of the random component of HRV, we utilized an autoregressive (AR) model (Akaike 1969) and measured the variance of residual time series reflecting the part of fluctuations that were not predicted by the AR model. To evaluate the age-dependent changes, we used the Allostatic State Mapping by Ambulatory ECG Repository (ALLSTAR) of 24-h electrocardiography (ECG) database, which consists of data collected from about 304 000 subjects. We randomly extracted data from subjects uniformly for all age groups who showed normal ECG and HRV profiles. Finally, as an additional pilot study, we investigated if the random component of HRV is low in patients with DD for their age.

This study was performed according to a protocol approved by the Ethics Review Committee of Nagoya City University Graduate School of Medical Sciences and Nagoya City University Hospital (No. 709 and No. 60160209).

Database

Clinical data

Data analysis

Preprocessing for HRV analysis

For both extracted and clinical 24-h ECG data, the temporal positions of all R-waves were detected, the rhythm annotations were given to all QRS complexes, and all errors in automated analysis were corrected manually by skilled clinical technologists. Only using R–R intervals consisting of two consecutive sinus rhythm beats (normal-to-normal [N–N] intervals), R–R time series data were then interpolated by a horizontal step function, i.e. the instantaneous R–R interval was defined as constant at the R–R interval value during the interval. Then, interpolated N–N interval time series were resampled for the entire 24 h and for nighttime data between 23:00 and 06:00 at 2 Hz for conventional HRV indices and at 1 Hz for AR model fitting. The nighttime data were used to investigate the possibility that age-dependent changes in HRV indices are the results of age-dependent changes in daily activities. If it applies to an HRV index, age-dependent changes in the index are expected to be lost or attenuated in nighttime data.

Assessment of conventional HRV indices

We calculated the conventional HRV measures recognized as mortality predictors in earlier studies (Camm et al 1996, Huikuri et al 2000, Bauer et al 2006, Kiyono et al 2008), which included SDNN, ultra-low-frequency (ULF, <0.0033 Hz), very-low-frequency (VLF, 0.0033–0.04 Hz), low-frequency (LF, 0.04–0.15 Hz), and high-frequency (HF, 0.15–0.40 Hz) component power, DC, short-term (α₁, 4–11 beats) and long-term (α₂, >11 beats) scaling exponents, spectral exponent β, and non-Gaussianity index λ_25s. The power of frequency components was calculated with a fast Fourier transform for the entire length of the N–N interval data as fluctuations in the respective frequency bands of individual components and was transformed into a natural logarithmic (ln) value. DC was measured as spontaneous increases within short (4 beats) overlapping segments of N–N intervals averaged over 24 h by the method of phase-rectified signal averaging (Kantelhardt et al 2007). Scaling exponents α₁ and α₂ were calculated by detrended fluctuation analysis (DFA) (Peng et al 1995) to quantify the correlation properties of short-term (4–11 beats) and long-term (>11 beats) fractal-like dynamics (Huikuri et al 2009). The spectral exponent β was calculated to examine 1/f^β power-law behavior of 24-h N–N interval variation according to the method reported previously (Bigger et al 1996, Huikuri et al 1998, Sakata et al 1999). Finally, λ_25s was calculated as an index of non-Gaussianity of the probability density function of the magnitude and occurrence frequency of abrupt intermittent changes in N–N intervals, and was interpreted as an index of sympathetic over-activity during daily life (Kiyono et al 2004, Kiyono et al 2008).

Assessment of random HRV component

To estimate the random component of HRV, we used an AR model fitting. By fitting an AR model, the fluctuation of time series can be separated into those of the regulated component (linear function of preceding values weighted by the AR coefficients) and the random component (residual of AR model), as shown in the equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x\left(t \right)=\sum\limits_{i=1}^{\,p}{a\left(i \right)x\left(t-i \right)+Z\left(t \right)},\nonumber \end{align}$$

where x(t) is the time series of N–N interval, a(i) is the AR coefficient, Z(t) is the time series of AR residual component, and p is the AR model order. Because Z(t) represents random fluctuations that cannot be predicted by the AR model, its variance $\sigma \left(Z\left(t \right) \right)$ can be considered as an estimate of the variance of the random component and $\sigma \left(x\left(t \right) \right)-\sigma (Z\left(t \right))$ as the estimate of the variance of the regulated component. We defined the ratio of variance of the random component in the total variance as AR variability (ARV):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle ARV=\frac{\sigma (Z\left(t \right))}{\sigma (x\left(t \right))}.\nonumber \end{align}$$

These methods were applied on the resampled N–N interval time series for both 24-h and nighttime (23:00–06:00) data. AR model fitting was performed by the algorithm published in Marple (1980). The model order was selected as that which minimized the final prediction error (FPE) figure of merit according to Akaike's information criterion (AIC) (Akaike 1969). Additionally, to examine the influence of the AR model order on ARV, ARV was calculated with model orders forcedly incremented from 2 to 300 by one in nine randomly selected subjects with different levels of ARV with the FPE criterion.

Statistical analysis

Sample extraction

Conventional HRV indices

Figure 1 shows age-dependent changes in conventional HRV indices. Heart rate was highest at the beginning of life, rapidly decreased from 0 to 20 years old, and gradually decreased after age 40. SDNN increased rapidly from 0 to 20 years old, decreased thereafter until age 40 and reached a plateau. DC increased from 0 to 15 years old in females and to 30 years old in males and decreased until age 80. Non-Gaussianity index λ_25s increased gradually after age 60. Power of frequency components showed a similar pattern, increasing from 0 to 15 years old and decreasing with advancing age. The power of LF and HF components was increased slightly at the end of life. Scaling exponent α₁ increased from 0 to 45 years old and decreased after age 75 in males, while it increased from 0 to 35 years old and decreased after age 65 in females. Scaling exponent α₂ decreased from 0 to 15 years old and increased with age until age 70 in males and age 80 in females. Finally, spectral exponent β increased from 10 years old until 70 years old and decreased thereafter in both sexes. Heart rate, HF, and α₂ were lower in males than in females (P  <  0.0001, P  =  0.004, and P  <  0.0001, respectively), while SDNN, DC, VLF, LF, LF/HF, α₁, and λ_25s were higher in males than in females (P  <  0.0001 for all but P  =  0.01 for SDNN).

Zoom In Zoom Out Reset image size

Random and regulated components of HRV

According to the FPE criterion, the median order (minimum, inter-quartile range, maximum) of the AR model for 24-h HRV was 189 (51, 142–267, 510) for males and 186 (46, 136–273, 501) for females. To examine the influence of the model order on estimated ARV, nine subjects with different ARV values were randomly selected and ARV was calculated with model orders forcedly incremented from 2 to 300 by 1. As shown in figure 2, while ARV decreased with increasing AR order, it saturated by the order of 40 and plateaued thereafter. Because the AR order with FPE criterion was  >45 in all subjects, we considered that the influence of the model order can be neglected.

Zoom In Zoom Out Reset image size

Figure 3 shows age-dependent changes in AR-model-derived HRV indices of 24-h and nighttime (23:00–06:00) data. Although nighttime data showed greater ARV and greater variance for the random component than those of 24-h data, the patterns of their age-dependent changes were similar, indicating that the age-dependent changes in 24-h data are not attributable to age-dependent changes in daily activities.

Zoom In Zoom Out Reset image size

The age-dependent changes in the regulated component were similar to those of SDNN in both sexes. On the other hand, the variance of the random component and ARV showed characteristic changes with age. The random component had a peak at 20 years old, decreased thereafter until 40, and increased again from age 75 in males and after age 85 in females. Furthermore, ARV was high at birth, decreased with age until 35 years old in men and 30 years old in women, and increased again from age 75 in men and age 85 in women. While the variance of the regulated component was greater in males (P  =  0.01), there was no significant gender difference in the variance of the random component or in ARV.

Table 1 shows the relationships between conventional and AR-model-derived HRV indices in 24-h data of both sexes (1930 males plus 1987 females). As expected, the variance of the regulated component was closely correlated with SDNN (r  =  0.96). ARV showed negative correlation with spectral exponent β (r  =  −0.57) and scaling exponent α₁ (r  =  −0.41) and mild positive correlations with HF (r  =  0.39).

Table 1. Correlation coefficients between conventional and AR-model-derived HRV indices in samples (1930 males and 1987 females) extracted from the ALLSTAR database.

	ARV	Random component	Regulated component
HR	0.07	−0.11	−0.30
SDNN	−0.17	0.28	0.96
DC	0.05	0.18	0.31
λ25s	−0.17	−0.13	0.11
ULF	−0.22	0.17	0.77
VLF	0.17	0.46	0.66
LF	0.33	0.50	0.52
HF	0.39	0.52	0.47
LF/HF	−0.20	−0.17	−0.07
Scaling exponent α1	−0.41	−0.31	0.02a
Scaling exponent α2	−0.30	−0.30	−0.06
Spectral exponent β	−0.57	−0.44	0.01a

^aExcept these, all other correlation coefficients are significant.

Patients with DD

In the present study, we found that the ratio of the random component of HRV estimated by ARV increases at the beginning and the end of life. Our findings are consistent with the hypothesis that the ratio of the random component of HRV increases during the developmental and degenerative stages of life, which may also support the concept of biological adaptation as attractor selection by the force of random fluctuation (Kashiwagi et al 2006).

The idea that an increase in random fluctuation associated with adaptation may be detected by AR model fitting has been deduced from the similarity of the AR model to the Langevin equation, which is used to represent the attractor selection model of biological adaptation (Kashiwagi et al 2006). A Langevin equation for the attractor selection model is written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{d}{dt}x=f\left(x \right)\cdot activity+\eta ,\nonumber \end{align}$$

where f(x) represents pulling force to the attractor, activity represents the force of bias that the system receives (goodness of state), and η represents the force of random fluctuation that the system receives. When organisms encounter a new environment or problem, they increase η to seek the new optimal attractors. By fitting an AR model,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x\left(t \right)=\sum\limits_{i=1}^{\,p}{a\left(i \right)x\left(t-i \right)+Z\left(t \right)},\nonumber \end{align}$$

the power of fluctuation can be separated into those of regulated and random fluctuation components. Because Z(t) represents random fluctuations that cannot be predicted by the model, its variance can be thought as an estimate of η. In a simulation study (not published), we generated simulated HRV signals by adding regular sinusoidal waves and random white noise in various ratios of their variances. By fitting an AR model to those simulated signals, ARV faithfully reflected the ratios of the variance of random noise to the total variance of the simulated signal.

In this study, we analyzed age-dependent changes in conventional HRV indices and compared them with those of the AR-model-derived indices. The age-dependent changes in SDNN, DC, and ULF, VLF and LF power, which increases during the early phase of life until 15–20 years old and declines thereafter, were close to those in the regulated component, indicating that most of the HRV consists of the regulated component. The age-dependent changes in HF power were somewhat similar to those of the random component and, consequently, there was mild correlation between them (r  =  0.52, respectively), suggesting associations between cardiac resting function (Hayano and Yasuma 2003) and random fluctuation of heart rate.

Interestingly, ARV showed the opposite age-dependent pattern to those of scaling exponent α₁ and spectral exponent β and negative correlations with them (r  =  −0.41 and  −0.57). An increase in α₁ above 1.0 is thought to represent a loss of complexity in short-term heart rate dynamics (Huikuri et al 2009). An increase in β towards 2 also represents a decrease in complexity and flexibility along with the changes from fractal dynamics to more correlated Brownian dynamics of fluctuation. These also indicate that 24-h HRV displaying dynamics closer to fractal noises yields a higher ratio of residual error in AR-model fitting, whereas 24-h HRV displaying dynamics closer to Brownian noises yields a lower ratio of residual error in AR fitting. Although these suggest that ARV may be related at least partially to the complexity and flexibility of heart rate dynamics, the pattern of age-dependent changes in ARV matched only modestly the patterns of these nonlinear HRV indices. In fact, the correlation coefficients indicate that age-dependent change in scaling and spectral exponents account for about 17% and 32%, respectively, of the age-dependent variations of ARV, suggesting that ARV may carry information that is not included in these nonlinear HRV indices.

We observed that analyses of 24-h and nighttime (23:00–06:00) data showed similar patterns of age-dependent changes in AR-model-derived HRV indices. Because data during sleep can be assumed to be less affected by individual differences in activities, the age-dependent changes observed in 24-h data do not seem attributable to the changes in daily activities with age. Furthermore, the variance of the random component was greater in nighttime data than in 24-h data despite the overall variance of R-R interval variability of nighttime data being about half of that in 24-h data. Consequently, ARV was  >1.5 times higher in nighttime than 24-h data. This indicates that the R-R interval variation was more random during nighttime, which seems consistent with earlier studies reporting increased nonlinearity (Beckers et al 2006) or complexity (Viola et al 2011) of heart rate dynamics during sleep. Given that the source of variations observed as HRV is derived through the autonomic nerves (Sands et al 1989), these findings indicate that the random behavior of heart rate dynamics is intrinsic to neural activities in the brain.

Finally, we observed that ARV was reduced in pediatric patients with DD for their age, although it should be considered as a preliminary finding of a pilot study. One may speculate, however, that the function of random fluctuation to facilitate attractor selections for development might be reduced in these patients. These patients also showed increased heart rate and scaling exponent α₁, suggesting increased tension during their daily life and less complexity of short-term heart rate dynamics.

This study has several limitations. First, the ECG data have been collected at a sampling frequency of 125 Hz. Although it is standard (125–128 Hz) in Holter ECG recordings for clinical uses, higher sampling rates might decrease the ratio of the random component due to smaller noises in the R–R interval measurement. Second, for all age groups from 0 to 100 years we used the same parameters of HRV analyses, such as the frequency bands for HRV components and the definitions of scaling regions for short- and long-term DFA. These frames themselves, however, may be changed with age and may have caused biases in the observed age-dependent changes in HRV indices. This possibility seems worth investigating in future studies. Third, the data in the ALLSTAR database has been collected for clinical reasons. In this sense, the ALLSTAR database may represent the characteristics of patients undergoing Holter ECG monitoring examination in Japan. Although we selected only subjects whose ECG showed normal sinus rhythm (>80% of all recorded beats), no abnormal finding of conduction or myocardial perfusion, and a low-risk profile from risk-predictive indices of HRV, we cannot exclude the possibility that the data might be affected by underlying diseases such as diabetes, hypertension, and cardiovascular diseases, and by possible differences in the prevalence of those diseases with age and sex. Nevertheless, the age-dependent changes we observed in conventional HRV indices are similar to those reported in an earlier study of 114 healthy subjects (Pikkujaemsae et al 1999), and the present study is first to demonstrate the complete profile of age- and sex-dependent changes in AR-model-derived HRV indices for the entire length of life. Finally, in this study we examined the hypothesis that the ratio of the random component of HRV increases during developmental and degenerative stages of life, based on the concept of biological adaptation as attractor selection by the force of random fluctuation. Our results supported the hypothesis, but at this time we cannot deny a logical leap in directly linking the results to the concept of biological adaptation. Our interpretation may be only one of several possibilities and should be positioned as one of the first challenges for a new interpretation of HRV.

Based on the concept that biological adaptation is regarded as attractor selection by the force of random fluctuation, we examined if the random component in HRV increases during developmental and degenerative stages of life. Our findings supported the hypothesis, indicating that the ratio of the random component of HRV is increased during the early and final phases of life and that it may be reduced in DD patients during their developmental stage.

This work was supported by a grant of the Knowledge Hub of Aichi, Japan (the Priority Research Project, P3-G1-S1-2b (J Hayano)) and the Japan Society for the Promotion of Science, Japan (Grant-in-Aid for Scientific Research (C) 23591055 (J Hayano); Grant-in-Aid for Scientific Research (C) 25461062 (H Fukuta); Grant-in-Aid for Scientific Research (B) 15H03095 (T Nakamura); and Grant-in-Aid for Scientific Research (A) 17H00878 (Y Yamamoto)).