Heart rate analysis in normal subjects of various age groups.

BACKGROUND: Analysis of heart rate variation (HRV) has become a popular noninvasive tool for assessing the activities of the autonomic nervous system (ANS). HRV analysis is based on the concept that fast fluctuations may specifically reflect changes of sympathetic and vagal activity. It shows that the structure generating the signal is not simply linear, but also involves nonlinear contributions. Linear parameters, Power spectral indice (LF/HF) is calculated with nonlinear indices Poincare plot geometry(SD1,SD2), Approximate Entropy (ApEn), Largest Lyapunov Exponent (LLE) and Detrended Fluctuation Analysis(DFA). The results show that, with aging the heart rate variability decreases. In this work, the ranges of all linear and nonlinear parameters for four age group normal subjects are presented with an accuracy of more than 89%. As a pre-analysis step, the HRV data is tested for nonlinearity using surrogate data analysis and the results exhibited a significant difference in the ApEn, LLE, SD1/SD2 and DFA parameters of the actual data and the surrogate data.
METHODS: The heart rate is analyzed using the various time domain parameters, frequency domain parameter and nonlinear parameters like Poincare geometry, ApEn, LLE and DFA.
RESULTS: In this work, the different linear and nonlinear parameters evaluated show a particular range for various cardiac abnormalities. And the results of these were subjected to 't' test with more than 89% confidence interval giving excellent 'p' values in all cases.
CONCLUSIONS: Heart rate variability (HRV) signal can be used as a reliable indicator of state of the heart. It becomes less random with the aging(less chaotic). This is evaluated by using various time domain, frequency domain and nonlinear parameters like SD1/SD2, ApEn, LLE alphas and alphal. Different ranges of non-linear parameters for various age groups are presented with 'p' value < or = 0.12.

Background

The Heart rate variability (HRV) is a non-invasive index of the neural control of the heart. HRV can be quantified by the simple calculation of the standard deviation of the R-R intervals. Furthermore, in the frequency domain, spectral analysis of HRV reveals three distinct frequency regions in the modulation of heart rate in humans. The typical spectral pattern, in normal conditions shows the presence of three frequency bands: a very low frequency (VLF) band from 0.00 to 0.03 Hz, a low frequency (LF) band from 0.03 to 0.15 Hz and a high frequency (HF) band in respiratory range generally more than 0.15 Hz. The power of LF component seems to be related to the vagal and sympathetic activities (LF component increases with every form of sympathetic stimulation), whereas the area of high frequency component (HF) provides a quantitative index of the influence of respiration on ECG signal and may be connected to the vagal activity. Thus LF/HF ratio is an important marker of sympathetic modulation or sympatho-vagal balance on heart rate variability control.

Physiological signals often vary in a complex and irregular manner. Analysis of linear statistics such as mean values, variability measures, and spectra of such signals generally does not address directly their complexity and thus may miss potentially useful information. Since the underlying mechanisms involved in the control of heart rate is mainly nonlinear, the application of nonlinear techniques seems appropriate [17,14,7,2].

Recently, new dynamic methods of HRV quantification have been used to uncover nonlinear fluctuations in heart rate that are not otherwise apparent. Several methods have been proposed: Lyapunov exponents [28], 1/f slope [16], approximate entropy(ApEn) [24] and detrended fluctuation analysis[22].

Compared to men, women are at lower risk of coronary heart disease [34] and of serious arrhythmias [3], suggesting a beneficial difference in autonomic control of heart rate.

Previous studies have assessed gender and age-related differences in time and frequency domain indices [27] and some nonlinear component of HRV [29]. There also seemed to be a significant difference between day and night hours when studying HRV indices using spectral and time domain methods [27,37].

Terry et al, have shown that the middle-aged women and men have a dominant parasympathetic and sympathetic regulation of the heart rate, respectively [31].

It is proved that, the heart rate variability depends on the sex also. The heart rate variability is more in the physically active young and old women [29,15]. The reduced HRV results in new cardiac events are studied by Hisako et al. It is proved by Emese et al, that the alert new borns heart rate variation is lower in the case of boys than in the case of girls [5]. The Heart rate variation for healthy subjects from 20–70 yrs is studied by Hendrik et al, and found that the HRV decreases with age and variation is more in the case of female than men [10]. Galeev et al, have analyzed the heart rate variation for healthy subjects of age from 6 to 16 yrs and observed the statistical and frequency domain variation with aging and gender [6].

In this work, we have analyzed the heart rate variation of normal subjects from 5–70 yrs in four groups (10 ± 5, 25 ± 5, 45 ± 5, 65 ± 5) by statistical method, frequency domain method and nonlinear method.

Materials and Method

In this project, ECGs are recorded from 150 healthy Chinese origin subjects (75 male and 75 female) with age ranging from 5 to 70 years for 20 minutes in lead II configuration. First the data is recorded in the relaxed sitting position. The total ECG data recorded is divided into four groups: little children, young, middle aged and old subjects. The number of subjects in each group is shown in table 1. Later for the same duration ECG is recorded with the same subject in the lying position. The ECG is recorded using the PowerLab/16SP system (ADI instrument). The data is sampled at a sampling rate of 400 Hz and stored in a random access file. The interval between two successive QRS complexes is defined as the r-r interval (tr-r seconds) [12] and the heart rate (beats per minute) is given as:

Table 1

Number of subjects in various groups

Age range	10 ± 5	25 ± 10	40 ± 15	60 ± 5
Number of subjects	25	50	40	35

HR = 60/tr-r     (1)

Table 1: Number of subjects in various groups

Analysis

The heart rate is analyzed in the time domain, frequency domain and by using nonlinear parameters.

Time-and Frequency-Domain Analysis of Heart Rate Variability

The time-and frequency-domain measures of HR variability were analyzed by the methods recommended by the Task Force of the European Society of Cardiology [30]. The standard deviation of all normal R-R intervals (SDNN) and the difference between the maximum hourly HR variability (circadian rhythm) were computed as standard time-domain measures of HR variability. Spectral power was quantified both by Fast Fourier Transform analysis and by autoregressive analysis in 4 frequency bands [11]:, ≤ 0.003 Hz(ultralow frequency [ULF]), 0.003 to 0.04 Hz(very-low frequency [VLF]), 0.04 to 0.15 Hz (low frequency [LF]), and 0.15 to 0.40 Hz (high frequency [HF]). ULF and VLF spectral components were computed over the entire recording interval by the fast Fourier method [11]. LF and HF components were computed from segments of 512 R-R intervals by the autoregressive method.

Poincare Plot Analysis

The Poincare plot, a technique taken from nonlinear dynamics, portrays the nature of R-R interval fluctuations. It is a plot in which each R-R interval is plotted as a function of the previous R-R interval. Poincare plot analysis is an emerging quantitative-visual technique whereby the shape of the plot is categorized into functional classes that indicate the degree of the heart failure in a subject [36] The plot provides summary information as well as detailed beat-to-beat information on the behavior of the heart [13].

The geometry of the Poincare plot is essential and can be described by fitting an ellipse to the graph. The ellipse is fitted onto the so called line-of-identity at 45° to the normal axis. The standard deviation of the points perpendicular to the line-of-identity denoted by SD1 describes short-term variability which is mainly caused by respiratory sinus arrhythmia (RSA). The standard deviation along the line-of-identity denoted by SD2 describes long-term variability.

Statistically, the plot displays the correlation between consecutive intervals in a graphical manner. Nonlinear dynamics considers the Poincare plot as the two dimensional (2-D) reconstructed R-R interval phase-spaces, which is a projection of the reconstructed attractor describing the dynamics of the cardiac system. The R-R interval Poincare plot typically appears as an elongated cloud of points oriented along the line-of-identity. The dispersion of points perpendicular to the line-of-identity reflects the level of short term variability. The dispersion of points along the line-of-identity is thought to indicate the level of long-term variability.

The Poincare plot may be analyzed quantitatively by calculating the standard deviations of the distances of the R-R(i) to the lines y = x and y = -x+2*R-Rm, where R-Rm is the mean of all R-R(i) [33]. The standard deviations are referred to as SD1 and SD2, respectively. SD1 related to the fast beat-to-beat variability in the data, while SD2 describes the longer-term variability of R-R(i) [33]. Figure 1 shows the Poincare plot of a normal young subject.

Figure 1

Poincare plot of a normal young subject.

Fig. 1 Poincare plot of a normal young subject

Approximate Entropy

Approximate entropy quantifies the regularity of time series. ApEn is also called a "regularity statistic". It is represented as a simple index for the overall "complexity" and "predictability" of each time series. In our study ApEn quantifies the regularity of the R-R interval. The more regular and predictable the R-R interval series, the lower will be the value of ApEn. On the other hand more the randomness in the R-R interval series, the higher will be the value of ApEn. Therefore, this method quantify the unpredictability of fluctuations in a time series such as an instantaneous R-R interval time series, R-R(i).

To compute ApEn of each data set, m-dimensional vector sequences [x[n]] were constructed from the R-R interval time series

x(n) = [R-R(n),.....R-R(N+m-1)]

Where the index n can take on values ranging from 1 to N-m+1 and N is the total number of data points in the R-R interval time series. If the distance between two vectors x(i) and x(j) is defined as d [x(i), x(j)], then we have

Where,

N, is a sequence of 1000 consecutive instantaneous R-R intervals.

m, specifies the pattern length which is 2 in this study.

r defines the criterion of similarity which has been set at 15% of the standard deviation of 1000 R-R intervals and (r) = { Number of x(j) such that d [x(i), x(j)] ≤ }/(N-m+1).

ApEn (N,r,m) is then defined as follows in equation (3)

ApEn = [ApEn(N,r,m) - ApEn(N,r,m + 1)]     (3)

The criterion of similarity, r, was chosen such that it was larger than most of the noise but at the same time not so large that detailed information about the system dynamics would be lost. On the basis of the work of Pincus et al [24-26], for m = 2, N = 1000 and r = 15% of the standard deviation of the data would produce reasonable statistical validity of ApEn.

In our study we use data set of 1000 adjacent R-R intervals. We divide the data set into smaller sets of length m = 2. This amounts to 500 smaller sub-sets. The next step is to determine the number of sub-sets that are within the criterion of similarity r = 15% of deviation of 1000 R-R intervals. Then we repeat the same process for the second sub-set till each sub-set is compared with the rest of the data set. This process computes the part of equation (2) and N-m+1 = 1000-2+1 = 999. We repeat the same process for m = 3. Approximate entropy is then calculated using equation (3).

Largest Lyapunov Exponent (LLE)

Lyapunov Exponent (λ) is a quantitative measure of the sensitive dependence on the initial conditions. It defines the average rate of divergence of two neighboring trajectories. An exponential divergence of initially nearby trajectories in phase space coupled with folding of trajectories, to ensure that the solutions will remain finite, is the general mechanism for generating deterministic randomness and unpredictability. Therefore, the existence of a positive λ for almost all initial conditions in a bounded dynamical system is widely used definition of deterministic chaos. To discriminate between chaotic dynamics and periodic signals Lyapunov Exponent (λ) are often used. It is a measure of the rate at which the trajectories separate one from other. The trajectories of chaotic signals in phase space follow typical patterns. Closely spaced trajectories converge and diverge exponentially, relative to each other. For dynamical systems, sensitivity to initial conditions is quantified by the Lyapunov Exponent (λ). They characterize the average rate of divergence of these neighboring trajectories. A negative exponent implies that the orbits approach a common fixed point. A zero exponent means the orbits maintain their relative positions; they are on a stable attractor. Finally, a positive exponent implies the orbits are on a chaotic attractor.

The algorithm proposed by Wolf et al [35] is used to Largest LE (LLE) from EEG data. For a given the time series x(t) for m dimensional phase space with delay coordinate t, that is a point on the attractor is given by

{ x(t), x(t + t),..., x(t + (m-1)t }

We locate nearest neighbor to initial point

{ x(t0), x(t0+ t),..., x(t0 + (m-1)t }

And denote the distance between these two points as L(t0). At a later time t1, initial length will evolve to length L'(t1). The mean exponential rate of divergence of two initially close orbits is characterized by

In implementation of this program, the following set of numerical parameters has to be chosen:

P = { m, t, T, Smax, Smin,thmax}

where m is the embedding dimension, t is delay, T being evaluation time (= t k+1 - t k-1) and Smax, Smin are the maximum and minimum separations of replacement point respectively and thmax is the maximum orientation error. According to Das et al [Das et al, 2002] an embedding dimension between 5 to 20 and a delay of 1 should be chosen when calculating LE for EEG data. In our analysis we have chosen an embedding dimension of 10 and delay of 1.

Detrended Fluctuation Analysis (DFA)

The concept of fractal is most often associated with geometrical objects satisfying two criteria: Self-similarity means that an object is composed of sub-units on multiple levels that statistically resemble the structure of the whole object. The second criterion for a fractal object is that it has a fractal dimension, also called fractal, that can be defined to be any curve or surface that is independent of scale. This concept of fractal structure can be extended to the analysis of physiological signals.

The DFA was used to quantify the fractal scaling properties of heart rate signals of short interval. Therefore, we apply DFA to the analysis of biological data, which calculates the root-mean-square fluctuation of integrated and detrended time series, permits the detection of intrinsic self-similarity embedded in a non-stationary time series, and also avoids the spurious detection of apparent self-similarity [1]. This method has been applied to a wide range of simulated and physiologic time series in recent years [19,8,9,21,22].

To illustrate the DFA algorithm, the total length of the heart rate signal (N) is integrated

Where HR(i) is the i th heart rate signal and HRis the average heart rate of N samples. Next, the integrated time series is divided into boxes of equal length n, a least-squares line is fit to the data (representing the trend in that box). The y-coordinate of the straight line segments is denoted by y(k). Then, we detrend the heart rate data y(k), by subtracting the local trend, y(k), in each box. The root mean square fluctuation of this integrated and detrended heart rate data is calculated by

This computation is repeated over all the box sizes to provide a relationship between F(n), the average fluctuation as a function of box size n.

In this study, the box size is ranged from 4 to ~300 beats. A box size larger than 300 beats would give a less accurate fluctuation value because of finite length effects of data.

Typically, F(n) will increase with box size. A linear relationship on a double-log graph indicates the presence of scaling i.e. F(n) ≈ nα. Under such conditions, the fluctuations can be characterized by a scaling exponent α, the slope of the line relating log F(n) to log(n). An α of 0.5 corresponds to white noise, α = 1 represents 1/f noise and α = 1.5 indicates Brownian noise or random walk. A good linear fit of the log F(n) to log(n) plot (DFA plot) indicates F(n) is proportional to nα, where α is the single exponent describing the correlation properties of the entire range of heart rate data [shown in Fig. 1]. However in some cases, we found that the DFA plot was not strictly linear but rather consisted of two distinct regions of different slopes separated at a break point n. This observation suggests there is a short range scaling exponent, α, over periods of 4 to 13 [n] beats, and a long – range exponent, α, over long periods [Woo et al, 1992; Kamen et al, 1996].

Results

Table 2 Variation of linear parameters for various age group normal subjects

Table 2

Variation of linear parameters for various age group normal subjects

PARAMETERS	10 ± 5 years	25 ± 10 years	40 ± 15 years	60 ± 5 years	'p' value
SDNN	92.96 ≤ 48.60357	81.30984 ≤ 38.87531	49.3888 ≤ 21.5635	70.89301 ≤ 61.20758	0.0007
SDSD	89.24 ≤ 64.86	80.56024 ≤ 62.02868	37.58513 ≤ 20.31865	78.50874 ≤ 93.46987	0.0067
RMSSD	88.482 ≤ 65.6591	80.1892 ≤ 62.338	37.56944 ≤ 20.31005	78.43348 ≤ 93.34064	0.008
Pnn50	12.9144 ≤ 8.469	15.2694 ≤ 9.97441	3.17521 ≤ 4.08388	9.77243 ≤ 19.8519	< 0.0001
Triangle index	0.018 ≤ 0.0049	0.0163 ≤ 0.0047	0.0123 ≤ 0.0046	0.017 ≤ 0.01482	0.0038
TiNN	121.51 ≤ 234	154.188 ≤ 372.983	421.60714 ≤ 619.8746193	113 ≤ 648	0.024

Table 3 Variation of parameter for various age group normal subjects in the frequency domain

Table 3

Variation of parameter for various age group normal subjects in the frequency domain

PARAMETERS	10 ± 5	25 ± 10	40 ± 15	60 ± 5	'p' value
LF/HF	1.425 ≤ 1.0591	1.26798 ≤ 0.88745	2.29766 ≤ 2.59557	1.57 ≤ 1.867	0.018

Table 4 Variation of nonlinear parameters for various age group normal subjects

Table 4

Variation of nonlinear parameters for various age group normal subjects

PARAMETERS	10 ± 5	25 ± 10	40 ± 15	60 ± 5	'p' value
SD1/SD2	0.5169 ≤ 0.193551	0.5265 ≤ 0.2357	0.4298 ≤ 0.191684	0.6055 ≤ 0.286834	0.12
ApEn	1.903 ≤ 0.34597	2.0722 ≤ 0.1849	1.8746 ≤ 0.3475	1.6837 ≤ 0.4083	< 0.0001
LLE	0.685206 ≤ 0.201294	0.6392875 ≤ 0.089508672	0.438793103 ≤ 0.177914783	0.434278 ≤ 0.294581	< 0.0001
αs	0.579 ≤ 0.35056	0.412118 ≤ 0.204044	0.138667 ≤ 0.216083	0.224479 ≤ 0.328036	< 0.0001
αl	2.0062 ≤ 0.33472	1.830926 ≤ 0.149596	1.674143 ≤ 0.205721	1.714868 ≤ 0.222235	< 0.0001

The results of the nonlinear parameters for four different age group is presented in table 2. To test for a statistical significance of difference in original SD1/SD2, ApEn, LLE, α, αand the surrogate data, surrogate data series were generated to match each original signal. Then it is subjected to Student t test distribution. We found that, the surrogate data for SD1/SD2, ApEn, LLE, α, αand their original data are distinct and they differ by more than 50%. This rejects the null hypothesis and hence the original data contain nonlinear features.

Surrogate Data

The purpose of surrogate data is to test for any nonlinearity in the original data [32]. To test if the attractor geometry and correlation dimension are truly due to chaotic dynamics one must examine these characteristics for surrogate data sets.

Nonlinear indexes such as ApEn are computed for several surrogate data series. Their values are compared with that assumed by the nonlinear index computed for the original index [32]. The demonstration of statistically significant difference in ApEn between the original and surrogate data are in keeping with the presence of nonlinear dynamics in the original data.

Surrogate data have Fourier decomposition with the same amplitudes as the empirical data decomposition but with random phase components. This is obtained from the Chaos Data Analyzer.

To test for a statistical significance of difference in original ApEn and the surrogate data, 10 surrogate data series were generated to match each original signal. Then it is subjected to Student t test distribution. We found that, the surrogate data ApEn and original data ApEn, are distinct and they differ by more than 50%. This rejects the null hypothesis and hence the original data contain nonlinear features.

Table 5 Results of various parameters using surrogate data

Table 5

Results of various parameters using surrogate data

PARAMETERS	Actual data	Surrogate Data	%Difference
SD1/SD2	0.5713	0.8652	51.4%
ApEn	1.4231	2.212	55%
LLE	0.4742	0.6379	34.5%
Alpha S	0.2134	0.3412	59.8%
Alpha L	1.722	2.345	36.1%

Discussion

Nikhil Iyengar et al, have applied DFA to the normal heart rate variability signal [18]. Recently, Echeverria et al, have applied detrended fluctuation analysis to the normal and Congestive Heart Failure (CHF) subjects. In this work, we have evaluated parameters α, α, ApEn and SD2 of the Poincare plots for the analysis of various cardiac abnormalities. And the results of these were subjected to't' test with more than 89% confidence interval giving excellent 'p' values in all cases.

The importance of ApEn lies in the fact that it is measure of the disorder in the heart rate signal. For the more regular and predictable heart rate values like complete heart block and ischemic/dilated cardiomyopathy type of abnormalities, the lower is the ApEn. On the other hand for the normal subjects where the heart rate is more random, ApEn has higher value. Hence the ApEn is used to find the unpredictability of fluctuations in the heart rate signals. This value decreases for the abnormal cardiac beats.

The pattern of Poincare plots of HRV data, its position and ranges of SD2 values are unique for particular type of cardiac abnormality. In this plot for the periodic data and for the normal heart rate signals the repeated stretching and folding of the trajectories, causes near by points to separate. Hence the position of the plot, its shape, SD2 computed with the HRV data may be considered as a characteristic parameter of diagnostic importance in clinical cardiology. For the normal subjects, shape of the plot is an ellipse and lies at the center of the quadrant and this position and shape of the plot shifts depending of the abnormality (Fig. 2-7). The value of SD2 indicates the long term variability. This value is incomparable with the SD1 in the case of abnormal cases.

Figure 2

F(n) plotted against several box sizes, n, on a log-log scale

The DFA technique is used to quantify the fractal scaling properties of short heart rate signals. This technique is a modification of root-mean-square analysis of random walks applied to non stationary data. From the root-mean-square fluctuation of an integrated and detrended time series is measured at different windows two parameters (αand α) derived. These parts have different range of values for various types of cardiac abnormalities (as indicated by very low 'p' value as shown in table 2). The short time correlation exists for abnormal signals with varying degrees and this correlation is less for normal subjects and for normal subject the long term correlation is more. In the normal case the short range correlation is less and long range correlation is more. As the signal becomes abnormal this short range correlation becomes more.

SD1/SD2 shows the ratio of short interval variation to the long interval variation. This ratio is more in the case of child and old male subjects, indicating lesser R-R variability in the small interval. This ratio is less in the middle aged subjects, indicating higher R-R variability.

The importance of ApEn lies in the fact that it is measure of the disorder in the heart rate signal. It is a measure, quantifying the regularity and complexity of time series. It has higher value in the case of children and young and this value decreases as the subject grows old. Hence, the ApEn will have smaller value for middle and old aged subjects, indicating smaller variability in the beat to beat.

Largest Lyapunov exponent's (LLE) quantify sensitivity of the system to initial conditions and gives a measure of predictability. This value decreases with the aging indicating that the heart rate variability becomes less chaotic as the healthy subject grows old.

Short range correlation αvalue is low in the abnormal signals. Similarly, the long range correlation αvalues are more for the normal signal indicating the fractalness in the data and this range decreases as the subject grows old.

Conclusion

Heart rate variability (HRV) signal can be used as a reliable indicator of state of the heart. It becomes less random with the aging(less chaotic). This is evaluated by using time domain, frequency domain and nonlinear parameters SD1/SD2, ApEn, LLE αand α. Different ranges of non-linear parameters for various age groups are presented with 'p' value ≤ 0.12.

Authors' Contributions

OWS and LYP carried out the data collection and software implementation. RAU, KN and CT carried out the design, analysis and coordination of the work. All authors read and approved the final manuscript.