Heart Rate Variability Analysis in Different Age and Pathological Conditions

Problem statement: Heart Rate Variability (HRV) has been used as a me sure of mortality primarily with patients who had undergone cardiac s urgery. The analysis of Heart Rate Variability (HRV) demands specific capabilities which are not p rovided either by parametric or nonparametric conventional estimation methods. The Empirical Mode Decomposition (EMD) adaptively estimates the Intrinsic Mode Functions (IMFs) of nonlinear nonsta tionary signals. Approach: The intrinsic mode functions estimated from the HRV signal were based on local characteristics of the signal. The princi ple objective was to analyze the HRV latencies of healt y subjects in different age and pathological conditions. The method was applied to HRV signal of 17 healthy young control subjects, 17 healthy old control subjects and 20 congestive heart failure pa tients for half hour duration. Results: The results showed that a healthy person’s HRV rapidly rises to its maximum response much earlier than the HRV of pathological subjects. The rising slope of the time scale’s plot discriminates the healthy controls an d pathological subjects with 100% sensitivity and spe cificity. Conclusion: This fact makes the method a promising approach to be applied in clinical practi ce as a screening test for specific risk-groups.


INTRODUCTION
Over the last 20 years there has been widespread interest in the study of variations in the beat-to-beat interval of heart known as heart rate variability or RR interval variations. Clinical depression strongly associated with mortality with such patients may be seen through a decrease in HRV. Heart rate is influenced by sympathetic and parasympathetic (vagal) activities of autonomous nervous system. The sympathetic activity accelerates the heart rate while the parasympathetic activity decelerates the heart rate. The influence of both branches of the autonomous nervous system is known as sympathovagal balance reflected in the HRV, which is a non invasive measure of the autonomous nervous system balance (Buccelletti et al., 2009;Feldman et al., 2010).
Empirical Mode Decomposition (EMD), introduced by Huang et al. (1998) is a method of decomposing nonlinear, non-stationary, multi component signals. The components resulting from EMD are called Intrinsic Mode Functions (IMFs). EMD is defined by an algorithm and has got no analytical formulation. Hence the decomposition is best understood by experimental investigation rather than analytical results. Being fully data dependent and highly adaptive it is found to be a highly efficient method of decomposing any nonlinear and nonstationary signals. Job Lindsen and Bhattacharya (2010) used EMD and Independent component analysis method to correct the blink artifacts (Lindsen and Bhattacharya, 2010). Ortiz et al. (2005) applied EMD method to decompose the fetal HRV series into its components in order to identify, the high frequency oscillations (Ortiz et al., 20005). Neto et al. (2004) applied EMD to situations where postural changes occur, provoking instantaneous changes in heart rate as a result of autonomic modifications. Shafqat et al. (2009) applied EMD to evaluate the effect of local anesthesia on HRV parameters. In this research the EMD method is used to analyze the HRV latencies of healthy subjects in different age and pathological conditions.

Empirical Mode Decomposition (EMD):
The Empirical mode decomposition is based on the following assumptions; 1.The signal has at least two extrema's-one maxima and one minima, 2. The characteristic time scale is defined as the time lapse between the extrema and 3. If the data is totally devoid of extrema, but has some inflection points, it can be differentiated one or more times to reveal the extrema. Final results are obtained by integration of the components. The essence of this method is to identify intrinsic oscillatory modes of their characteristic time scales in the data empirically and then decompose the data accordingly. The components obtained as a result of the decomposition process are termed as Intrinsic Mode Functions (IMFs). Formally an IMF can be defined as a function that satisfies the following two conditions 1.In the whole data set, the number of extrema and the number of zero crossings must be either equal to or differ by at most one and 2.At any point the mean value of the envelope defined by local minima and the envelope defined by local maxima is zero. The name IMF is adopted since it represents the oscillation mode imbedded in the data. With this definition of IMF in each cycle, defined by the zero crossings involves only one mode of oscillation with no complex riding waves.

Datasets used in the analysis:
To study the intrinsic mode functions of HRV in different age and pathological condition, half an hour duration HRV signal from three different groups of subjects were considered for the analysis: • 17 healthy young control subjects without any clinical evidence of heart disease • 17 healthy old control subjects without any clinical evidence of heart disease • 20 Congestive Heart Failure (CHF) The ECG data for the three groups has been collected from the biomedical website, http://www.physionet.org. The healthy subjects ECG data was drawn from the Fantasia database and the CHF data from the BIDMC-CHFDB.
Discrete event series, R i -R i-1 intervals as a function of R i occurrence times, was constructed by an adaptive QRS detector algorithm. The QRS detector was based on the one presented Christov (2004). As a result of the detection algorithm, an unevenly sampled RR interval time series was obtained. In order to recover an evenly sampled signal from the irregularly sampled event series, cubic interpolation was applied.

Methodology:
The EMD method is a sifting process that estimates the local time scales of HRV signal. It involves the following steps, leading to a decomposition of the signal S(t) into its constituents components: • x (an auxillary variable) is set to the signal S(t) to be analyzed and a variable k, which is the number of estimated IMFs, is set to zero • Splines are fitted to the upper extrema and the lower extrema. This will define the lower (LE) and Upper Envelopes (UE) • The average envelope, m, is calculated as the arithmetic mean between UE and LE • A candidate IMF, h, is estimated as the difference between x and m • If h does not fulfill the criteria defining an IMF, it is assigned to the variable x and the steps (b)-(d) are repeated. Otherwise, if h is an IMF then the procedure moves to step (f) • If h is an IMF it is saved as c k , where k is the kth component • The mean squared error, mse, between two consecutive IMFs c k-i and c k is calculated and this value is compared to a stopping condition (usually a very small value, i.e. 10 -5 ) • If the stopping condition is not reached, the partial residue, r k , is estimated as the difference between a previous partial residue r k-1 and c k and its content is assigned to the dummy variable x and the steps of (b)-(d) are repeated • If the stopping condition is reached then the sifting process is finalized and the final residue r final can be estimated as the difference between S(t) and the sum of all IMFs The criterion used to state whether h is an IMF or not is to verify whether h satisfies the two conditions that define an IMF. Currently, there is no set of equations to estimate IMFs; therefore, the sifting procedure described above is an empirical technique, is employed for this purpose. An example of typical IMF is shown in Fig. 1.
When the sifting process stops, the original signal S(t) can be represented as:

Fig. 1: A typical IMF
Where n is the number of IMFs, c k the kth IMF and r final is the final residue.
Equation 1 indicates that IMFs can be linearly combined in order to obtain the decomposed signal S(t).

Time domain HRV measures:
The power of the nth IMF was computed as: where, c n is the nth IMF and j=1,….N samples and average period (mean period) of the IMF, c n is: where, dist refers to the distance between the first and last zero crossings and Zc is the number of zero crossings.
Interpretation of intrinsic mode function: By the nature of the decomposition procedure, the signal is decomposed into N fundamental components, each with distinct time scale. More specifically, the first component is the smallest time scale component which corresponds to the highest frequency component of the data. As the decomposition process proceeds, the time scale increases and hence, the mean frequency of the scales decreases. The average period of intrinsic time scale almost doubles that of the previous one, suggesting that the EMD behaves like a dyadic filter. Based on this a time scale filtering method can be build up as: When l=1 and h < N, it is a high-pass filter; when l > 1 and h=N, it is a low-pass filter; when 1<l≤ h< N, it is a band-pass filter.

RESULTS
The method is applied to half an hour duration HRV measurements of 17 healthy young controls, 17 healthy old controls and 20 congestive heart failure patients. At first, the method is demonstrated by applying it to a typical young and old Control and a Heart Failure (CHF) subject. Afterwards, the discrimination of the three subject groups with respect to different parameters obtained from the method will be shown. The typical signals, IMFs and the reconstructed signals are shown in Fig. 2. For the typical young and old control subject and CHF patient, records f1y07 and f1o04 and CHF 02 records from physionet website were selected, shown in (Fig. 2 a(i),b(i) andc(i)). The method adaptively decomposes the three signals into IMFs as shown in Fig. 2. a(ii),b(ii) andc(ii). The average periods, absolute powers and normalized powers of IMFs are computed. The signal is reconstructed using the IMFs leaving out the residue component were shown in (Fig.2.a (iii), b(iii) andc(iii)).
The absolute powers (V n ) and average periods (T n ) of the IMFs are computed using simple formulae given in equation (2) and (3). The average period of the successive IMFs almost doubles the previous IMF's average period. The residue part of the signal is a monotonic trend with average period zero (T n =0).
The computed average periods (T n ) of IMFs for the 3 subjects are given in Table 1.
The method decomposes the healthy young and old control's HRV signal (f1y07 and f1o04) into eight IMFs and CHF patient's HRV into ten IMFs. The additional component in CHF patient's HRV was due to the latencies present in the signal. Plotting the average periods of IMFs against its IMF number gives an exponential graph as shown in Fig. 3. The average period of IMFs of CHF 02 subject was significantly lower (Table 1) in value and the rate of increase (slope) w.r.to IMFs also smaller compared to healthy controls. From Fig. 3, the rising slopes of the exponential curves of healthy controls are significantly higher compared to CHF patients' curve. The slope of the exponential curve was approximated using y=ab x , a curve fitting mathematical equation. The variable 'y' represents the average period of IMFs and 'x' represents the IMF numbers. Variable 'a' represents the slope of the exponential graph and 'b' is the scaling factor (approximately 2) of the time scales. The parameters 'a' and 'b' are estimated using simple least square curve fitting technique (Ramana,0000). The estimated 'a' and 'b' parameters for the typical healthy young and old controls and CHF patient are given in Table 2.  The slope 'a' of the average period curve of CHF patient is significantly lower (0.99042) compared to healthy controls (2.1683 and 1.6244). The average periods of successive IMFs of healthy controls and CHF patients approximately doubles (1.94 to 2.07 times) as that of the previous IMFs average period, as defined by the parameter 'b'.
The computed absolute power (V n ) of IMFs for the 3 subjects were presented in Fig. 4. For healthy young control subject (f1y07) the absolute power was high in all IMFs. Two significant peaks are observed at IMF 1 and IMF 5. For healthy old control subject (f1o04) the power is less compared to healthy young in all IMFs except in IMF4 and dominates all the other IMFs. But for CHF 02 the absolute powers of all IMFs were completely suppressed.
The IMFs powers can be calculated in normalized units (n.u) which represent the relative value of each IMF power in proportion to the total power. Normalization tends to minimize the effect on the powers of IMFs, if the total power changes. The computed normalized powers (Nn) of IMFs for the 3 subjects were presented in Fig. 5 for comparison. In healthy young control subject (f1y07) the relative powers of all IMFs (except IMF 7and8) are high. In healthy old control subject (f1o04) the relative power of IMF 4 is very dominating compared to all IMFs. But for CHF 02 the relative powers of all IMFs were suppressed except for the last (IMFs 9 and 10) two IMFs power, the last scale (IMF 10) is dominating compared to all the IMFs. Suppression of power in the lower scales (higher frequencies) makes the system less adaptive.
The contribution of first six IMFs of healthy young and old subjects to the HRV measurements are approximately 92 % and 94%. The contribution of third, fourth and fifth IMFs of healthy old is approximately 72 %. But, the contribution of first six IMFs of CHF patient to the total power was approximately 50%. The results show that more power in the lowest (IMF 1) scale increases the power of high frequency component in healthy young subject. More power in the middle (IMF 4, 5) scales increases the power of lower frequency components in healthy old subjects. In CHF subjects the power in the lowest scale and middle scales (high and lower frequency components) are much suppressed, but the power in the last scales (lowest frequency components) were higher. The latencies of the IMFs can be presented visually by plotting the cumulative sum of IMF's normalized powers against it's IMF number as shown in Fig. 6.
The curve for healthy young control subject is more rapidly rising from the first IMF and approaches the maximum earlier. The curve for healthy old control subjects is initially very slow in the first three IMFs but after the third IMF more rapidly rising and approaches the maximum earlier. But, the curve for CHF subject is slow rising and approaches the maximum later. These observations are very characteristics for all the three groups of signals.
The group's average plots of average periods are shown in Fig. 7. The group's average plots show that the slopes of the healthy controls are higher than CHF group patients. The slopes, parameter 'a' completely discriminates the healthy controls and CHF group patient with 100% sensitivity and specificity for a threshold value of 1.4 as given in Fig. 11. The estimated parameters are given in Table-3 and 4. The scaling parameter 'b' is approximately equal to two as given in Fig. 12.     The group's absolute power average plot (Fig. 8) shows the absolute power of young subjects are significantly high compared to old and CHF group patients. In healthy young subjects the power in the lowest scale (highest frequency) dominates other scales.   In healthy old subjects the power in the middle scale (IMF 4) dominates other scales. In CHF group the power is completely suppressed in all scales. But in the CHF group's normalized power average plot (Fig. 9) the normalized powers in the last scales (IMFs 9 and10) are very dominating compared to all lower scales. The cumulative sum of normalized powers of IMFs average plot (Fig. 10) completely discriminates the healthy controls and patient groups in the middle scales (4 th , 5 th and 6 th IMFs). The CHF patients average curve is significantly very slow in rising and approaches the maximum later. The estimated parameters 'a and b' of exponential curves for all the 54 records are given in the following Tables (3 and 4). Parameter 'a' completely discriminates (as in Fig.11) the healthy and patient's groups and the plots of parameter 'b' shows the method is almost like a dyadic filter (as in Fig.12).

DISCUSSION
In this study, a practical method for analyzing the HRV latencies is presented. It was shown that the latencies of HRV signal discriminates the healthy subjects and congestive heart failure subjects significantly. As a specific application the method was applied to half an hour HRV measurement of healthy controls and congestive heart failure patients and a good discrimination of the two groups was obtained.
The EMD method estimates the local time scales adaptively which reflects the intrinsic properties of the signal. The first six scales (1-6 IMFs) of healthy control subjects contribute significant power (more than 85% of the total power) to the signal total power, but for congestive heart failure system the power in these scales are much suppressed (less than 50% of total power) as shown in the Fig. 13. This feature makes the healthy systems to reach its maximum response much earlier and makes the system more adaptive than congestive heart failure patients.

CONCLUSION
The common hypothesis is that the human cardiovascular system is a highly complex adaptive system and that the complexity of its behavior allows for the broadest range of adaptive responses. The proposed technique is simple and adaptive method to analyze the complex HRV signal. The fastness in reaching maximum response of the healthy system represents its more adaptiveness for particular level of input and the slowness in reaching maximum response (more latency) of CHF subjects represents the system's inability to respond quickly for various levels of inputs. All patients and healthy controls are discriminated correctly with 100% sensitivity and specificity. This fact makes the method a promising approach to be applied in clinical practice as a screening test for specific risk-groups.