- Email: [email protected]
Investigating fractal property and respiratory modulation of human heartbeat time series using empirical mode decomposition
Medical Engineering & Physics 32 (2010) 490–496
Contents lists available at ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Investigating fractal property and respiratory modulation of human heartbeat time series using empirical mode decomposition Jia-Rong Yeh a , Wei-Zen Sun b , Jiann-Shing Shieh a,∗ , Norden E. Huang c a b c
Department of Mechanical Engineering, Yuan Ze University, 135 Yuan-Tung Rd., Chung-Li, Taoyuan, Taiwan Department of Anaesthesiology, College of Medicine, National Taiwan University, Taipei, Taiwan Research Center for Adaptive Data Analysis, National Central University, Taoyuan, Taiwan
a r t i c l e
i n f o
Article history: Received 25 September 2009 Received in revised form 24 February 2010 Accepted 26 February 2010 Keywords: Heartbeat interval Long-term correlation Empirical mode decomposition Intrinsic mode function Intrinsic mode analysis Respiratory sinus arrhythmia
a b s t r a c t The human heartbeat interval reflects a complicated composition with different underlying modulations and the reactions against environmental inputs. As a result, the human heartbeat interval is a complex time series and its complexity can be scaled using various physical quantifications, such as the property of long-term correlation in detrended fluctuation analysis (DFA). Recently, empirical mode decomposition (EMD) has been shown to be a dyadic filter bank resembling those involved in wavelet decomposition. Moreover, the hierarchy of the extracted modes may be exploited for getting access to the Hurst exponent, which also reflects the property of long-term correlation for a stochastic time series. In this paper, we present significant findings for the dynamic properties of human heartbeat time series by EMD. According to our results, EMD provides a more accurate access to long-term correlation than Hurst exponent does. Moreover, the first intrinsic mode function (IMF 1) is an indicator of orderliness, which reflects the modulation of respiratory sinus arrhythmia (RSA) for healthy subjects or performs a characteristic component similar to that decomposed from a stochastic time series for subjects with congestive heart failure (CHF) and atrial fibrillation (AF). In addition, the averaged amplitude of IMF 1 acts as a parameter of RSA modulation, which reflects significantly negative correlation with aging. These findings lead us to a better understanding of the cardiac system. © 2010 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction It has been established that the heartbeat interval is a complex time series, and so quantifying the complexity of such signals could be a powerful tool to understand the underlying controlling mechanisms of human heartbeat and other related situations of physiological systems. Previous investigations on the dynamic properties of human heartbeat time series have been conducted using various approaches, such as the entropy in multi-scale entropy (MSE) [1] and the fractal property in detrended fluctuation analysis (DFA) [2,3]. These physical properties were used to quantify the complexity of human heartbeat with some degree of success. Recently, the fractal complexity measure ‘˛’ of DFA has further been proven to have clear diagnostic and prognostic value [3]. The DFA of heart rate variability was presented as a time-domain technique, in which the time series of RR intervals is cumulatively summed and then cut into short segments. Within each segment, the degree to which the cumulative time series is dispersed away
∗ Corresponding author. Tel.: +886 3 4638800x2470; fax: +886 3 4558013. E-mail address: [email protected] (J.-S. Shieh).
from its linear trend is measured [4]. However, the calculation of DFA is based on a linear analysis technique, in which the trends of different timescales are derived by the linear method of least square line. In addition, since the fractal property exists only in a limited range of timescales, the calculation of DFA is used only to quantify fractal properties of the integrated time series, but not the original time series. On contrast to the linearly detrending technique of the least squares line used in DFA, empirical mode decomposition (EMD) has been proposed as an adaptive algorithm for decomposing intrinsic mode functions (IMF) from nonlinear signals [5]. This also functions as a detrending technique that can be used to obtain intrinsic trends adaptive to the nature of signals [6]. Recently, a combination of EMD and fractal dimension (FD) analysis was proposed to form a denoising EMD-FD filter for diagnosing explosive lung sounds [7]. Moreover, further studies also proved that EMD acts essentially as a dyadic filter bank [8,9], and the hierarchy of the extracted mode functions may be similarly exploited for getting access to the Hurst exponent [9]. Therefore, we here suggest a new scaling factor to quantify the fractal property of human heartbeat time series in this EMD-based analysis algorithm. In addition, we also found an extra parameter, different from the fractal property of human heartbeat time series in our previous
1350-4533/$ – see front matter © 2010 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2010.02.022
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
approach using EMD-based analysis [10]. In the current approach, the timescale index performs as an extra parameter, which reflects the distribution of intrinsic timescales of IMFs decomposed from human heartbeat time series by EMD. It shows the major differences between heartbeat time series of healthy subjects and subjects with heart disease. Moreover, since EMD functions as a dyadic filter bank for a time series with fractal properties, the sequence of intrinsic timescales of IMFs decomposed from human heartbeat time series is a geometric series with a ratio of around 2. Therefore, to determine the fundamental modulation of human cardiac system can help understand the underlying mechanism that dominates the distribution of intrinsic timescales. As shown in related studies, respiratory sinus arrhythmia (RSA) modulation is a dominant mechanism of cardiac system [11–13], and RSA can be easily observed in the heartbeat time series measured from healthy subjects. Moreover, EMD appropriately acts as a powerful tool to extract the particular component of a human heartbeat time series, that reflects the basic rhythm of the cardiac system. The present study has two primary goals. The first is to present a new scaling exponent for the fractal property of a human heartbeat time series using EMD. For the purpose to examine the performance of the EMD-based scaling exponent of fractal property, the analysis results using a new scaling exponent were compared with those by rescaled range analysis. However, since the scale of fractal property for human heartbeat time series was not a confirmable value for a real-world signal, it is hard to evaluate the accuracy of scaling factor for fractal property of human heartbeat time series. Therefore, simulated time series of Gaussian noise with different fractal properties were used to examine the accuracy of proposed new scaling factor. Pearson’s correlation coefficient was used to present the consistency between the values from the scaling exponent and the original fractal properties for the simulated time series of Gaussian noise. Moreover, the analysis result of the fractal property of human heartbeat time series from EMD-based analysis was compared with those from the previous analysis algorithm of DFA. Secondly, we wish to detect the respiratory modulation of the human cardiac system by an advanced analysis algorithm of EMD. Integrating the analysis results of fractal property and RSA modulation of human heartbeat time series can provide a better understanding for the underlying control mechanisms of the human cardiac system. 2. Material In this study, a database downloaded from Physiobank [14] is used as the study material. This database includes 40 healthy subjects with subgroups of young and elderly (20 young and 20 elderly), 43 subjects with severe congestive heart failure (CHF), and 9 subjects with atrial fibrillation (AF). The ages of subjects are from 21 to 34 (25.95 ± 4.31) years old for the subgroup of young and from 68 to 85 (74.55 ± 4.45) years old for the subgroup of elderly. This database has been used in many previous studies of human heartbeat analysis using different analysis algorithms. It provides samples of human heartbeat time series for four groups (i.e., healthy young, healthy elderly, subjects with CHF and subjects with AF). Therefore, the analysis results of this proposed EMD-based algorithm can be compared with results of previous studies to illustrate its advantages. 3. Empirical mode decomposition (EMD) Empirical mode decomposition (EMD) uses the mean of envelopes, which reflect the highest level of local trends, to derive the local oscillation riding on the local trend. A decomposed local oscillation is called an intrinsic mode function (IMF). The signal is considered to consist of all levels of local oscillations and the final
491
monotonic trend. A residue is considered to be a trend that consists of all un-decomposed local oscillations (i.e., IMFs) and the final monotonic trend. For a signal x(t), the proposed algorithm of EMD is suggested as follows [5]: 1) Connect the sequential local maxima (respective minima) to derive the upper (respective lower) envelop using cubic spline. 2) Derive the mean of the envelope, m(t), by averaging the upper and lower envelopes. 3) Extract the temporary local oscillation h(t) = x(t) − m(t). 4) Repeat steps (1)–(3) on the temporary local oscillation h(t) until m(t) is close to zero. Then, h(t) is an IMF noted as c(t). 5) Compute the residue r(t) = x(t) − c(t). 6) Repeat the steps (1) through (5) using r(t) for x(t) to generate the next IMF and residue. Therefore, the original signal x(t) can be reconstructed by the following formula: x(t) =
n
ci (t) + rn (t)
(1)
i=1
where ci (t) is the ith IMF (i.e., local oscillation), n is the number of IMFs, and rn (t) is the nth residue (i.e., local trend). 4. Estimating fractal property of nonlinear time series using EMD-based analysis 4.1. Averaged period and energy density used in the period-energy plot According to previous studies, EMD acts essentially as a dyadic filter bank to decompose a stochastic time series in wavelet decompositions [15]. Each residual from decomposition is the intrinsic trend with its corresponding intrinsic timescale (i.e., averaged period of the decomposed IMF). Therefore, the plot of energy densities against their corresponding averaged periods for IMFs decomposed by EMD reflects the power-law correlation between energy and timescale. The slope of the energy–timescale plot acts as a scaling factor of the fractal property. Here, energy density and averaged period of each IMF can be calculated by the following equations: 1 2 [cn (j)] N N
En =
(2)
j=1
T¯ n =
Sln T,n d ln T
(Sln T,n /T )d ln T
(3)
where En is the energy density of the nth IMF, N is the number of samples, n is the number of IMFs, Sln T,n is the Fourier spectrum of the nth IMF as a function of ln T, T is the period, and T¯ n is the averaged period of the nth IMF. For the purpose to demonstrate the calculating process of averaged period and energy density, the heartbeat time series of an AF patient was decomposed into the first 5 IMFs by EMD to derive the Fourier spectra and the energy–timescale plot as shown in Fig. 1. Fig. 1(a) shows the original heartbeat time series; Fig. 1(b) shows the first 5 IMFs decomposed from the original heartbeat time series; Fig. 1(c) shows the Fourier spectra of IMFs 1–5; Fig. 1(d) shows the energy–timescale plot using energy densities and averaged periods of IMFs 1–5. Furthermore, we defined a new scaling exponent of the fractal property using the slope of the energy–timescale plot in EMDbased analysis. As is well known, white noise is a fractal Gaussian
492
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
Fig. 1. Illustration of the calculating process of energy–timescale plot using the heartbeat time series of an AF patient: (a) the original heartbeat time series; (b) the first 5 IMFs decomposed from the original heartbeat time series; (c) the Fourier spectra of IMF 1–5; (d) the logarithmic energy–timescale plot using energy densities (ln E) and intrinsic timescale (ln T).
noise with Hurst exponent H of 0.5. 0.5 < H < 1.0 reflects positive long-term correlation and 0 < H < 0.5 reflects a negative long-term correlation [16,17]. Moreover, the Wood and Chan algorithm [18] has previously been applied to generate a simulated time series of fractal Gaussian noise with H varying from 0.1 to 0.9. Therefore, we conducted two numerical experiments to explore the accuracy of fractal property estimation using EMD-based analysis. In this first numerical experiment, 1000 simulated time series with N = 4096 for each fractal Gaussian noise with H varying from 0.1 to 0.9 using step 0.2 were generated by the Wood and Chan algorithm. The first 5 IMFs of the simulated time series of fractal Gaussian noises were decomposed to show the characteristics of fractal Gaussian noises derived by EMD. The ensemble results of EMD-based analysis are shown in Fig. 2 using averages and standard deviations. Slopes of the energy–timescale plots for fractal Gaussian noises with various H perform as a good indicator for fractal property of Gaussian noise. There, we defined the scaling exponent of fractal property in the presented EMD-based analysis. In order to derive a positive value for a new scale in EMD-based analysis for a fractal Gaussian noise with positive long-term correlation, the new scale was defined as the slope of energy–timescale plot plus 1 (since the slope of energy–timescale plot is −1 for white noise). In addition, we had conducted a second numerical experiment using 1000 simulated time series (N = 4096) of fractal Gaussian noise with random H on the range between 0.1 and 0.9. Then both rescaled range analysis (R/S analysis) [19] and EMD-based analysis were used to estimate the fractal properties of those simulated time series. Fig. 3 shows the distributions of the new scale vs. the original H of simulated time series and values of Hurst exponent evaluated by R/S analysis vs. the original H. Pearson’s correlation coefficient [20] was used estimate the accuracy of new scale in EMD-based analysis.
4.2. Comparison between the estimations of fractal property using DFA and EMD-based analysis As mentioned above, the new scale was defined to quantify the fractal property of nonlinear time series in the EMD-based analysis. The fractal property is an important characteristic of nonlinear signals since it reflects the functional mechanism of auto-correlation for nonlinear systems. In previous studies, detrended fluctuation analysis (DFA) has been used to investigate the fractal property
Fig. 2. Overall results of energy–timescale plots derived by the first 5 IMFs in the EMD-based analysis for fGns with various Hurst exponents from 0.1 to 0.9, where log T is the logarithmic intrinsic timescale and log E is the logarithmic energy density.
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
493
Fig. 3. Scaling exponents by rescaled range analysis and scaling factor of EMD-based analysis using the slope of energy–timescale plot derived from the first 5 IMFs for the simulated time series of fGn with varying Hurst exponent: (a) scaling exponent of rescaled range analysis; (b) scaling factor of EMD-based analysis.
of human heartbeat time series within short-term and long-term time scales [3,4]. The short-term scaling exponent (˛1 ) within a timescale range from 3 to 11 is a good physiological index for vagal modulation and the autonomous nerve system (ANS) for subjects with healthy functions of cardiac systems. It fits for ANS activity monitoring with some limitations. In contrast to short-term autocorrelation, the overall scaling exponent (˛) of DFA, which reflects the physiological mechanism of long-term auto-correlation, has been applied in functional verification of human cardiac system in the pioneering approach using DFA [2]. This study aimed to verify a human cardiac system evaluation using the long-term fractal property of heartbeat time series. Therefore, we compared analysis results for the long-term fractal property of human heartbeat time series by both DFA-based and EMD-based analyses. Fig. 4 shows the statistical analysis results of four groups (i.e., healthy young, healthy elderly, CHF and AF) in mean and standard deviation derived by DFA-based and EMD-based analyses. Here, Pearson’s correlation coefficient was applied to quantify the correlation between the analysis results by DFA-based and EMD-based analyses. The high value of Pearson’s correlation coefficient (r = 0.799) shows high consistency between results from both analysis algorithms of long-term auto-correlation. Statistical analysis by one way analysis of variance (ANOVA), P-value < 0.001 shows significant differences among long-term fractal properties of human heartbeat time series for the four groups. But long-term correlation is not the only property in human heartbeat analysis by EMD-based algorithm. And, according to our previous study, there is another property dependent on the distribution of intrinsic timescales in EMD-based analysis [10].
5. Underlying modulation of respiratory sinus arrhythmia (RSA) As is well known, modulation of RSA is an important auto-regulatory mechanism of the human cardiac systems. The modulation of RSA can be observed easily in a controlledrespiration experiment. In this experiment, a young healthy volunteer breathed following the instruction with a fixed rhythm (10 s per cycle) given by an mp3 voice recorder. Fig. 5 shows the heartbeat time series in two different situations (i.e., controlled and un-controlled breathing). In the situation of controlled breath-
ing, RR intervals significantly fluctuated with the same rhythm of breathing because the influence of RSA modulation had been magnified by controlled breathing. In contrast to data from controlled breathing, RR intervals were also affected by the RSA modulation but it is more complicated when the subject breaths naturally. Therefore, RSA is a significant and basic modulation of the cardiac system and it should be an intrinsic component of human heartbeat time series.
5.1. Classification of subjects with or without heart disease using RSA modulation According to our previous approach to human heartbeat analysis using EMD-based DFA, the distribution of intrinsic timescales derived by EMD is a good indicator for the underlying physiological situation of the cardiac system, verifying whether or not there is heart disease. Moreover, EMD also acts as a dyadic filter bank for human heartbeat time series. Therefore, the sequence of intrinsic timescales derived by EMD is a geometric series with ratio of 2 that is dependent on the first intrinsic timescale (i.e., averaged period of IMF 1). In order to find the intrinsic component of a heartbeat time series that reflects RSA modulation of the cardiac system, we calculated the intrinsic timescales (i.e., averaged periods of IMFs) for four groups, as shown in Table 1. We found that averaged periods of IMF 1 for two healthy groups (i.e., young and elderly) are similar to the rhythm of human respiration, which is in the range from 12 to 20 cycles per minute (equal to 3–5 s per cycle). But the intrinsic timescale of IMF 1 for the two groups with heart diseases (i.e., CHF and AF) is different from the rhythm of human respiration. They have an intrinsic timescale similar to IMF 1 decomposed from a totally stochastic signal, which has an averaged period of 2.81 sampling intervals (i.e., heartbeat interval) of the time series. For human heartbeat time series with heartbeat intervals from 0.5 to 1.0 s (i.e., heart rate from 60 to 120 beats per minute at rest), the time period of 2.81 heartbeat intervals is in the range from 1.405 to 2.81 s. Moreover, most subjects with CHF or AF have a high heart rate (average rate of around 100 beats per minute at rest). For a totally stochastic heartbeat time series with averaged heart rate of 100 beats per minute, the average period of IMF 1 is 1.686 s. Then, we found that the intrinsic timescale of IMF 1 for a healthy sub-
494
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
Fig. 4. Analysis results of new scaling factor of EMD-based analysis and DFA scaling exponent. Pearson’s correlation coefficient between results using two different parameters is 0.799: (a) new scaling factor of EMD-based analysis; (b) DFA scaling exponent.
Fig. 5. Human heartbeat time series in breath-control experiment.
ject has a rhythm similar to respiration (between 3 and 5 s) but the timescale for a subject with CHF or AF has rhythm similar to IMF 1 decomposed from a totally stochastic time series (close to 1.68 s). Finally, according to the results shown in Table 1, the statistical difference between intrinsic timescales of healthy group (i.e., young or elderly) and group with heart disease (i.e., CHF or AF) is significant (P-value < 10−4 by t-test). But, the statistical difference between groups young and elderly is insignificant (P-value > 0.05 by t-test). Also, the statistical difference between groups CHF and AF is insignificant (P-value > 0.05 by t-test). Therefore, an intrinsic timescale of IMF 1 decomposed from a human heartbeat time series can be a clue to evaluate the RSA modulation of a cardiac system. In a healthy subject, the RSA modulation is clear and reflected by the intrinsic timescale of IMF 1. But, the RSA modulation is diminished and approaches to the characteristic
component decomposed from a stochastic time series for a subject with CHF or AF. Based on this finding, an intrinsic timescale of IMF 1 can be used as an indicator for the orderliness of a human heartbeat time series. A low value of intrinsic timescale similar to IMF 1 decomposed from a stochastic time series reflects a high degree of disorderliness, as characteristic of a subject with CHF or AF, whereas a high value indicates RSA modulation of a healthy subject. In addition, the receiver operating characteristic (ROC) curve can be used to examine the sensitivity and specificity of classification using different parameters. Therefore, the ROC curve was used to estimate the sensitivity and specificity of classification of subjects with or without heart disease. This used the both fractal property and intrinsic timescale of IMF 1 to show the critical parameter for classification. Fig. 6 shows that the intrinsic timescale of IMF 1 is a better indicator for distinguishing healthy subjects from subjects
Table 1 Statistical values for intrinsic timescale of IMF 1 for heartbeat intervals of four groups (i.e., healthy young, healthy elderly, CHF and AF). Group
Number of cases
Age
Young Elderly CHF AF
20 20 43 9
25.95 74.55 55.37 43.48
± ± ± ±
4.31 4.45 11.32 12.25
Severity
Intrinsic timescale of IMF 1 (s)
N/A N/A NYHA classes I, II, & III N/A
3.34 3.28 1.98 1.77
± ± ± ±
0.45 0.55 0.31 0.12
Notation Similar to the rhythm of respiration Similar to the rhythm of respiration Similar to IMF 1 decomposed from a random time series Similar to IMF 1 decomposed from a random time series
N/A means not available and NYHA is the abbreviation of New York Heart Association. Here, NYHA means the classification of New York Heart Association.
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
495
inverse relationship with aging showing that averaged amplitude of IMF 1 has a negative connection with aging. However, since the aging process is exceedingly complex, many variables should be taken into account in the experimental protocol for a further study, which aims to understand the underlying mechanisms of aging. Fig. 7 presents the statistical results of intrinsic timescale and averaged amplitude for four groups. Here, intrinsic timescales of healthy groups (i.e., young and elderly) are significantly different from those of groups with heart disease and averaged amplitudes of groups of healthy young and AF are significantly different from those of groups of elderly and CHF. However, IMF 1 reflects different underlying modulations for different groups. For healthy groups, it reflects the strength of RSA modulation. For example, the cardiac system of a healthy young subject has a stronger RSA modulation than that of a healthy elderly subject. The averaged amplitude of IMF 1 for a young subject is larger than that for an elderly subject. Therefore, the average amplitude of IMF 1 has a positive correlation with the strength of RSA modulation and a negative correlation with aging. Fig. 6. ROC curves for classification of subjects with or without heart disease using two different parameters of the fractal property (i.e., new scaling factor of EMDbased analysis) and averaged period of IMF 1.
with CHF or AF. As the physiological meaning, RSA modulation is an important characteristic of human cardiac system. Thus, intrinsic timescale of IMF 1 is an effective indicator in diagnosis of heart disease (i.e., CHF and AF in this study). 5.2. Connection between aging and RSA modulation As mentioned above, IMF 1 of human heartbeat time series for a healthy subject reflects RSA modulation of the cardiac system. An intrinsic timescale of IMF 1 was used to check with the rhythm of respiration as an indicator for subjects with or without heart disease. Furthermore, to understand the connection between aging and RSA modulation, we conducted a comparison between averaged amplitudes of IMF 1 for groups of healthy young and healthy elderly. These two groups have a significantly statistical difference in average amplitude of IMF 1 (P-value = 0.0039 by t-test), with the average amplitude of IMF 1 for the young group significantly bigger than that for the elderly group. Thus, RSA modulation has an
6. Discussion and conclusions In this study, we present a new approach for human heartbeat analysis, using an EMD-based analysis algorithm. We extended a finding based on the function of EMD, in which the hierarchy of the extracted mode functions may be exploited for obtaining the fractal property of time series, to develop a new scaling parameter to quantify the long-term fractal property of human heartbeat time series. According to the results of a numerical experiment, the new parameter had been proven to perform as a more accurate scaling index than the Hurst exponent by R/S analysis via the simulated time series of fractal Gaussian noise with various Hurst exponents from 0.1 to 0.9. Moreover, its results are highly consistent with those using the long-term scaling exponent by DFA in human heartbeat analysis. Theoretically, fractal property reflects the characteristic of selfsimilarity for a signal among continuous observing scales. For a stochastic time series, both DFA and the EMD-based analysis reflect a consistent characteristic of self-similarity by straight lines of logarithmic energy/timescale plots. However, a human heartbeat time series is not a totally stochastic time series, which contains dominant components caused by underling modulations. Hence, the
Fig. 7. Statistical values of intrinsic timescales and averaged amplitudes of IMF 1 for four groups.
496
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
logarithmic period-energy plots of EMD-based analysis are not a well-defined straight line for heartbeat time series of healthy subjects and some CHF patients. On contrast to EMD-based analysis, the trends are approximated by least square lines in DFA, which are insignificantly affected by internal modulation. Therefore, the analysis results of EMD-based analysis are more variable than those of DFA. Excluding concerns of the internal modulation in human heartbeat time series, the fractal property of human heartbeat time series is an important aspect for system identification of human cardiac system, consistent with the results of previous studies of human heartbeat analysis. Fractal property performs well to show the characteristic property of a human’s heartbeat time series for verifying an ill subject with CHF or AF, but it does not provide a critical parameter in system verification between subjects with or without heart disease. In addition, the statistical difference between the fractal properties of groups of young and elderly is insignificants. On the other hand, previous studies also showed that shortterm scaling exponent (˛1 ) of DFA works to quantify the ensemble responses of RSA and ANS modulations [21,22]. However, the application of short-term scaling exponent of DFA in human heartbeat analysis has some limitations. A basic limitation is that ˛1 works only for healthy subjects (i.e., subjects without heart disease). Thus, EMD had been used to investigate the intrinsic modulation in human heartbeat time series. According to our findings in this EMDbased analysis, IMF 1 reflects the RSA modulation of cardiac systems only for healthy subjects but not for subjects with CHF or AF. For subjects with CHF or AF, IMF 1 has similarly characteristic timescale to that of IMF 1 decomposed from a stochastic time series. This shows that the human heartbeat time series of subjects with CHF or AF contains high-frequency components, which diminishes the component of RSA modulation. This causes that short-term scaling exponent of DFA works well only for quantifying the RSA and ANS modulation of healthy subjects but not for subjects with CHF or AF. Furthermore, according to our new findings, the intrinsic timescale of IMF 1 can be used as an index of orderliness for human heartbeat time series. A value of intrinsic timescale of IMF 1 similar to the rhythm of respiration indicates the underlying orderliness of RSA modulation to the human cardiac system. On contrast to the orderliness caused by RSA modulation, IMF 1 with an intrinsic timescale of 1.68 s presents a high degree of randomness for human heartbeat time series. This finding is the most important contribution of our study. Simply speaking, RSA modulation is a key factor for diagnosis of subjects with or without heart disease (i.e., CHF or AF in this study) and the intrinsic timescale of IMF 1 decomposed from human heartbeat time series is a key indicator for RSA modulation. Moreover, average amplitude of IMF 1 for subjects with heart disease reflects the effects caused by the presence of noise. The present noise in human heartbeat time could be met in the high-frequency part and embedded in IMF 1. The noise present in heartbeat time series is comparatively weak to result small amplitude of IMF 1 for a CHF patient. But, the present noise dominates the heartbeat time series and produces large amplitude of IMF 1 for an AF patient. Therefore, amplitude of IMF 1 reflects the strength of noise presenting in heartbeat time series as a key parameter for distinguishing the subject with CHF or AF. In summary, our previous study showed there is another dimension, different from the fractal property evaluated by the original scaling exponent of DFA, which is a new parameter for the distribution of intrinsic timescales. In this study, we found that RSA contributes to the fundamental modulation of the human cardiac system in healthy subjects. Thus, irregularity of RSA modulation,
as shown in IMF 1, is a significant sign to a cardiac system with CHF or AF. Furthermore, EMD is a powerful tool to extract the RSA modulation from human heartbeat time series, which is the most important contribution in the EMD-based analysis. Acknowledgements The authors wish to thank Prof. Ary L. Goldberger and Dr. C.K. Peng (the director and co-director of the Rey Institute for Nonlinear Dynamics in Medicine at the Beth Israel Deaconess Medical Center of Harvard Medical School) for valuable discussion. We gratefully acknowledge support from the National Science Council (NSC) of Taiwan (Grant Number NSC96-2221-E-155-015-MY3-2) for supporting this research. Conflict of interest statement The authors want to disclose there are no potential conflicts of interest with other people or organizations in this study. References [1] Costa M, Goldberger AL, Peng CK. Multiscale entropy analysis of complex physiologic time series. Phys Rev Lett 2002;89(6):068102. [2] Goldberger AL, Peng CK, Lipsitz LA. What is physiologic complexity and how does it change with aging and disease? Neurobiol Aging 2002;23:23–6. [3] Ho KK, Moody GB, Peng CK, Mietus JE, Larson MG, Levy D, et al. Predicting survival in heart failure case and control subjects by use of fully automated methods for deriving nonlinear and conventional indices of heart rate dynamics. Circulation 1997;96:842–8. [4] Francis DP, Willson K, Georgiadou P, Wensel R, Davies LC, Coats A, et al. Quantitative general theory for periodic breathing in chronic heart failure and its clinical implications. Circulation 2000;102:2214–21. [5] Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc Roy Soc Lond A 1998;454:903–95. [6] Wu Z, Huang NE, Long ST, Peng CK. On the trend, detrending, and variability of nonlinear and nonstationary time series. PNAS 2007;104(38):14889–94. [7] Hadjileontiadis LJ. A novel technique for denoising explosive lung sounds empirical mode decomposition and fractal dimension filter. IEEE Eng Med Biol Mag 2007;26(1):30–9. [8] Wu Z, Huang NE. A study of the characteristics of white noise using the empirical mode decomposition method. Proc Roy Soc Lond A 2004;460:1597–611. [9] Flandrin P, Rilling G, Goncalces P. Empirical mode decomposition as a filter bank. IEEE Signal Process Lett 2004;11:112–4. [10] Yeh JR, Fan SZ, Shieh JS. Human heart beat analysis using a modified algorithm of detrended fluctuation analysis based on empirical mode decomposition. Med Eng Phys 2009;31:92–100. [11] Hirsch JA, Bishop B. Respiratory sinus arrhythmia in humans: how breathing pattern modulates heart rate. Am J Physiol Heart Circ Physio 1981;241:H620–9. [12] Berntson GG, Carcioppo JT, Quigley KS. Cardiac psychophysiology and autonomic pace in humans: empirical perspectives and conceptual implications. Psychol Bull 1993;114(2):296–322. [13] Craft N, Schwartz JB. Effects of age on intrinsic heart rate, heart rate variability, and AV conduction in healthy humans. Am J Physiol 1995;268:H1441–52. [14] Databases are available at, http://www.physionet.org/, See Goldberger AL, et al. Circulation 2000;101:E215. [15] Flandrin P, Goncalves P. Empirical mode decomposition as a data-driven wavelet-like expansions. Int J Wavelet Multires Info Proc 2004;2:477–96. [16] Mandelbrot BB, van Ness JW. Fractional Brownian motions, fractional noises and applications. SIAM Rev 1968;10:422–37. [17] Hurst HE. Long-term storage capacity of reservoirs. Trans Am Soc Civil Eng 1951;116:770. [18] Wood AT, Chan G. Simulation of stationary processes in [0,1]d . J Comput Graph Stat 1994;3:409–32. [19] Alessio E, Carbone A, Castelli G, Frappiettro V. Second-order moving average and scaling of stochastic time series. Eur Phys J B 2002;27:197–200. [20] Glantz SA. Primer of biostatistics. 6th ed. Singapore: McGraw-Hill Inc.; 2005. [21] Peng CK, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationary heart-beat time series. Chaos 1995;5:82–7. [22] Mourot L, Boudaddi M, Gandelin E, Cappelle S, Nguyen NU, Wolf JP, et al. Conditions of autonomic reciprocal interplay versus autonomic co-activation: effects on non-linear heart rate dynamics. Autonom Neurosci: Basic Clin 2007;137:27–36.
Contents lists available at ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Investigating fractal property and respiratory modulation of human heartbeat time series using empirical mode decomposition Jia-Rong Yeh a , Wei-Zen Sun b , Jiann-Shing Shieh a,∗ , Norden E. Huang c a b c
Department of Mechanical Engineering, Yuan Ze University, 135 Yuan-Tung Rd., Chung-Li, Taoyuan, Taiwan Department of Anaesthesiology, College of Medicine, National Taiwan University, Taipei, Taiwan Research Center for Adaptive Data Analysis, National Central University, Taoyuan, Taiwan
a r t i c l e
i n f o
Article history: Received 25 September 2009 Received in revised form 24 February 2010 Accepted 26 February 2010 Keywords: Heartbeat interval Long-term correlation Empirical mode decomposition Intrinsic mode function Intrinsic mode analysis Respiratory sinus arrhythmia
a b s t r a c t The human heartbeat interval reflects a complicated composition with different underlying modulations and the reactions against environmental inputs. As a result, the human heartbeat interval is a complex time series and its complexity can be scaled using various physical quantifications, such as the property of long-term correlation in detrended fluctuation analysis (DFA). Recently, empirical mode decomposition (EMD) has been shown to be a dyadic filter bank resembling those involved in wavelet decomposition. Moreover, the hierarchy of the extracted modes may be exploited for getting access to the Hurst exponent, which also reflects the property of long-term correlation for a stochastic time series. In this paper, we present significant findings for the dynamic properties of human heartbeat time series by EMD. According to our results, EMD provides a more accurate access to long-term correlation than Hurst exponent does. Moreover, the first intrinsic mode function (IMF 1) is an indicator of orderliness, which reflects the modulation of respiratory sinus arrhythmia (RSA) for healthy subjects or performs a characteristic component similar to that decomposed from a stochastic time series for subjects with congestive heart failure (CHF) and atrial fibrillation (AF). In addition, the averaged amplitude of IMF 1 acts as a parameter of RSA modulation, which reflects significantly negative correlation with aging. These findings lead us to a better understanding of the cardiac system. © 2010 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction It has been established that the heartbeat interval is a complex time series, and so quantifying the complexity of such signals could be a powerful tool to understand the underlying controlling mechanisms of human heartbeat and other related situations of physiological systems. Previous investigations on the dynamic properties of human heartbeat time series have been conducted using various approaches, such as the entropy in multi-scale entropy (MSE) [1] and the fractal property in detrended fluctuation analysis (DFA) [2,3]. These physical properties were used to quantify the complexity of human heartbeat with some degree of success. Recently, the fractal complexity measure ‘˛’ of DFA has further been proven to have clear diagnostic and prognostic value [3]. The DFA of heart rate variability was presented as a time-domain technique, in which the time series of RR intervals is cumulatively summed and then cut into short segments. Within each segment, the degree to which the cumulative time series is dispersed away
∗ Corresponding author. Tel.: +886 3 4638800x2470; fax: +886 3 4558013. E-mail address: [email protected] (J.-S. Shieh).
from its linear trend is measured [4]. However, the calculation of DFA is based on a linear analysis technique, in which the trends of different timescales are derived by the linear method of least square line. In addition, since the fractal property exists only in a limited range of timescales, the calculation of DFA is used only to quantify fractal properties of the integrated time series, but not the original time series. On contrast to the linearly detrending technique of the least squares line used in DFA, empirical mode decomposition (EMD) has been proposed as an adaptive algorithm for decomposing intrinsic mode functions (IMF) from nonlinear signals [5]. This also functions as a detrending technique that can be used to obtain intrinsic trends adaptive to the nature of signals [6]. Recently, a combination of EMD and fractal dimension (FD) analysis was proposed to form a denoising EMD-FD filter for diagnosing explosive lung sounds [7]. Moreover, further studies also proved that EMD acts essentially as a dyadic filter bank [8,9], and the hierarchy of the extracted mode functions may be similarly exploited for getting access to the Hurst exponent [9]. Therefore, we here suggest a new scaling factor to quantify the fractal property of human heartbeat time series in this EMD-based analysis algorithm. In addition, we also found an extra parameter, different from the fractal property of human heartbeat time series in our previous
1350-4533/$ – see front matter © 2010 IPEM. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.medengphy.2010.02.022
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
approach using EMD-based analysis [10]. In the current approach, the timescale index performs as an extra parameter, which reflects the distribution of intrinsic timescales of IMFs decomposed from human heartbeat time series by EMD. It shows the major differences between heartbeat time series of healthy subjects and subjects with heart disease. Moreover, since EMD functions as a dyadic filter bank for a time series with fractal properties, the sequence of intrinsic timescales of IMFs decomposed from human heartbeat time series is a geometric series with a ratio of around 2. Therefore, to determine the fundamental modulation of human cardiac system can help understand the underlying mechanism that dominates the distribution of intrinsic timescales. As shown in related studies, respiratory sinus arrhythmia (RSA) modulation is a dominant mechanism of cardiac system [11–13], and RSA can be easily observed in the heartbeat time series measured from healthy subjects. Moreover, EMD appropriately acts as a powerful tool to extract the particular component of a human heartbeat time series, that reflects the basic rhythm of the cardiac system. The present study has two primary goals. The first is to present a new scaling exponent for the fractal property of a human heartbeat time series using EMD. For the purpose to examine the performance of the EMD-based scaling exponent of fractal property, the analysis results using a new scaling exponent were compared with those by rescaled range analysis. However, since the scale of fractal property for human heartbeat time series was not a confirmable value for a real-world signal, it is hard to evaluate the accuracy of scaling factor for fractal property of human heartbeat time series. Therefore, simulated time series of Gaussian noise with different fractal properties were used to examine the accuracy of proposed new scaling factor. Pearson’s correlation coefficient was used to present the consistency between the values from the scaling exponent and the original fractal properties for the simulated time series of Gaussian noise. Moreover, the analysis result of the fractal property of human heartbeat time series from EMD-based analysis was compared with those from the previous analysis algorithm of DFA. Secondly, we wish to detect the respiratory modulation of the human cardiac system by an advanced analysis algorithm of EMD. Integrating the analysis results of fractal property and RSA modulation of human heartbeat time series can provide a better understanding for the underlying control mechanisms of the human cardiac system. 2. Material In this study, a database downloaded from Physiobank [14] is used as the study material. This database includes 40 healthy subjects with subgroups of young and elderly (20 young and 20 elderly), 43 subjects with severe congestive heart failure (CHF), and 9 subjects with atrial fibrillation (AF). The ages of subjects are from 21 to 34 (25.95 ± 4.31) years old for the subgroup of young and from 68 to 85 (74.55 ± 4.45) years old for the subgroup of elderly. This database has been used in many previous studies of human heartbeat analysis using different analysis algorithms. It provides samples of human heartbeat time series for four groups (i.e., healthy young, healthy elderly, subjects with CHF and subjects with AF). Therefore, the analysis results of this proposed EMD-based algorithm can be compared with results of previous studies to illustrate its advantages. 3. Empirical mode decomposition (EMD) Empirical mode decomposition (EMD) uses the mean of envelopes, which reflect the highest level of local trends, to derive the local oscillation riding on the local trend. A decomposed local oscillation is called an intrinsic mode function (IMF). The signal is considered to consist of all levels of local oscillations and the final
491
monotonic trend. A residue is considered to be a trend that consists of all un-decomposed local oscillations (i.e., IMFs) and the final monotonic trend. For a signal x(t), the proposed algorithm of EMD is suggested as follows [5]: 1) Connect the sequential local maxima (respective minima) to derive the upper (respective lower) envelop using cubic spline. 2) Derive the mean of the envelope, m(t), by averaging the upper and lower envelopes. 3) Extract the temporary local oscillation h(t) = x(t) − m(t). 4) Repeat steps (1)–(3) on the temporary local oscillation h(t) until m(t) is close to zero. Then, h(t) is an IMF noted as c(t). 5) Compute the residue r(t) = x(t) − c(t). 6) Repeat the steps (1) through (5) using r(t) for x(t) to generate the next IMF and residue. Therefore, the original signal x(t) can be reconstructed by the following formula: x(t) =
n
ci (t) + rn (t)
(1)
i=1
where ci (t) is the ith IMF (i.e., local oscillation), n is the number of IMFs, and rn (t) is the nth residue (i.e., local trend). 4. Estimating fractal property of nonlinear time series using EMD-based analysis 4.1. Averaged period and energy density used in the period-energy plot According to previous studies, EMD acts essentially as a dyadic filter bank to decompose a stochastic time series in wavelet decompositions [15]. Each residual from decomposition is the intrinsic trend with its corresponding intrinsic timescale (i.e., averaged period of the decomposed IMF). Therefore, the plot of energy densities against their corresponding averaged periods for IMFs decomposed by EMD reflects the power-law correlation between energy and timescale. The slope of the energy–timescale plot acts as a scaling factor of the fractal property. Here, energy density and averaged period of each IMF can be calculated by the following equations: 1 2 [cn (j)] N N
En =
(2)
j=1
T¯ n =
Sln T,n d ln T
(Sln T,n /T )d ln T
(3)
where En is the energy density of the nth IMF, N is the number of samples, n is the number of IMFs, Sln T,n is the Fourier spectrum of the nth IMF as a function of ln T, T is the period, and T¯ n is the averaged period of the nth IMF. For the purpose to demonstrate the calculating process of averaged period and energy density, the heartbeat time series of an AF patient was decomposed into the first 5 IMFs by EMD to derive the Fourier spectra and the energy–timescale plot as shown in Fig. 1. Fig. 1(a) shows the original heartbeat time series; Fig. 1(b) shows the first 5 IMFs decomposed from the original heartbeat time series; Fig. 1(c) shows the Fourier spectra of IMFs 1–5; Fig. 1(d) shows the energy–timescale plot using energy densities and averaged periods of IMFs 1–5. Furthermore, we defined a new scaling exponent of the fractal property using the slope of the energy–timescale plot in EMDbased analysis. As is well known, white noise is a fractal Gaussian
492
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
Fig. 1. Illustration of the calculating process of energy–timescale plot using the heartbeat time series of an AF patient: (a) the original heartbeat time series; (b) the first 5 IMFs decomposed from the original heartbeat time series; (c) the Fourier spectra of IMF 1–5; (d) the logarithmic energy–timescale plot using energy densities (ln E) and intrinsic timescale (ln T).
noise with Hurst exponent H of 0.5. 0.5 < H < 1.0 reflects positive long-term correlation and 0 < H < 0.5 reflects a negative long-term correlation [16,17]. Moreover, the Wood and Chan algorithm [18] has previously been applied to generate a simulated time series of fractal Gaussian noise with H varying from 0.1 to 0.9. Therefore, we conducted two numerical experiments to explore the accuracy of fractal property estimation using EMD-based analysis. In this first numerical experiment, 1000 simulated time series with N = 4096 for each fractal Gaussian noise with H varying from 0.1 to 0.9 using step 0.2 were generated by the Wood and Chan algorithm. The first 5 IMFs of the simulated time series of fractal Gaussian noises were decomposed to show the characteristics of fractal Gaussian noises derived by EMD. The ensemble results of EMD-based analysis are shown in Fig. 2 using averages and standard deviations. Slopes of the energy–timescale plots for fractal Gaussian noises with various H perform as a good indicator for fractal property of Gaussian noise. There, we defined the scaling exponent of fractal property in the presented EMD-based analysis. In order to derive a positive value for a new scale in EMD-based analysis for a fractal Gaussian noise with positive long-term correlation, the new scale was defined as the slope of energy–timescale plot plus 1 (since the slope of energy–timescale plot is −1 for white noise). In addition, we had conducted a second numerical experiment using 1000 simulated time series (N = 4096) of fractal Gaussian noise with random H on the range between 0.1 and 0.9. Then both rescaled range analysis (R/S analysis) [19] and EMD-based analysis were used to estimate the fractal properties of those simulated time series. Fig. 3 shows the distributions of the new scale vs. the original H of simulated time series and values of Hurst exponent evaluated by R/S analysis vs. the original H. Pearson’s correlation coefficient [20] was used estimate the accuracy of new scale in EMD-based analysis.
4.2. Comparison between the estimations of fractal property using DFA and EMD-based analysis As mentioned above, the new scale was defined to quantify the fractal property of nonlinear time series in the EMD-based analysis. The fractal property is an important characteristic of nonlinear signals since it reflects the functional mechanism of auto-correlation for nonlinear systems. In previous studies, detrended fluctuation analysis (DFA) has been used to investigate the fractal property
Fig. 2. Overall results of energy–timescale plots derived by the first 5 IMFs in the EMD-based analysis for fGns with various Hurst exponents from 0.1 to 0.9, where log T is the logarithmic intrinsic timescale and log E is the logarithmic energy density.
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
493
Fig. 3. Scaling exponents by rescaled range analysis and scaling factor of EMD-based analysis using the slope of energy–timescale plot derived from the first 5 IMFs for the simulated time series of fGn with varying Hurst exponent: (a) scaling exponent of rescaled range analysis; (b) scaling factor of EMD-based analysis.
of human heartbeat time series within short-term and long-term time scales [3,4]. The short-term scaling exponent (˛1 ) within a timescale range from 3 to 11 is a good physiological index for vagal modulation and the autonomous nerve system (ANS) for subjects with healthy functions of cardiac systems. It fits for ANS activity monitoring with some limitations. In contrast to short-term autocorrelation, the overall scaling exponent (˛) of DFA, which reflects the physiological mechanism of long-term auto-correlation, has been applied in functional verification of human cardiac system in the pioneering approach using DFA [2]. This study aimed to verify a human cardiac system evaluation using the long-term fractal property of heartbeat time series. Therefore, we compared analysis results for the long-term fractal property of human heartbeat time series by both DFA-based and EMD-based analyses. Fig. 4 shows the statistical analysis results of four groups (i.e., healthy young, healthy elderly, CHF and AF) in mean and standard deviation derived by DFA-based and EMD-based analyses. Here, Pearson’s correlation coefficient was applied to quantify the correlation between the analysis results by DFA-based and EMD-based analyses. The high value of Pearson’s correlation coefficient (r = 0.799) shows high consistency between results from both analysis algorithms of long-term auto-correlation. Statistical analysis by one way analysis of variance (ANOVA), P-value < 0.001 shows significant differences among long-term fractal properties of human heartbeat time series for the four groups. But long-term correlation is not the only property in human heartbeat analysis by EMD-based algorithm. And, according to our previous study, there is another property dependent on the distribution of intrinsic timescales in EMD-based analysis [10].
5. Underlying modulation of respiratory sinus arrhythmia (RSA) As is well known, modulation of RSA is an important auto-regulatory mechanism of the human cardiac systems. The modulation of RSA can be observed easily in a controlledrespiration experiment. In this experiment, a young healthy volunteer breathed following the instruction with a fixed rhythm (10 s per cycle) given by an mp3 voice recorder. Fig. 5 shows the heartbeat time series in two different situations (i.e., controlled and un-controlled breathing). In the situation of controlled breath-
ing, RR intervals significantly fluctuated with the same rhythm of breathing because the influence of RSA modulation had been magnified by controlled breathing. In contrast to data from controlled breathing, RR intervals were also affected by the RSA modulation but it is more complicated when the subject breaths naturally. Therefore, RSA is a significant and basic modulation of the cardiac system and it should be an intrinsic component of human heartbeat time series.
5.1. Classification of subjects with or without heart disease using RSA modulation According to our previous approach to human heartbeat analysis using EMD-based DFA, the distribution of intrinsic timescales derived by EMD is a good indicator for the underlying physiological situation of the cardiac system, verifying whether or not there is heart disease. Moreover, EMD also acts as a dyadic filter bank for human heartbeat time series. Therefore, the sequence of intrinsic timescales derived by EMD is a geometric series with ratio of 2 that is dependent on the first intrinsic timescale (i.e., averaged period of IMF 1). In order to find the intrinsic component of a heartbeat time series that reflects RSA modulation of the cardiac system, we calculated the intrinsic timescales (i.e., averaged periods of IMFs) for four groups, as shown in Table 1. We found that averaged periods of IMF 1 for two healthy groups (i.e., young and elderly) are similar to the rhythm of human respiration, which is in the range from 12 to 20 cycles per minute (equal to 3–5 s per cycle). But the intrinsic timescale of IMF 1 for the two groups with heart diseases (i.e., CHF and AF) is different from the rhythm of human respiration. They have an intrinsic timescale similar to IMF 1 decomposed from a totally stochastic signal, which has an averaged period of 2.81 sampling intervals (i.e., heartbeat interval) of the time series. For human heartbeat time series with heartbeat intervals from 0.5 to 1.0 s (i.e., heart rate from 60 to 120 beats per minute at rest), the time period of 2.81 heartbeat intervals is in the range from 1.405 to 2.81 s. Moreover, most subjects with CHF or AF have a high heart rate (average rate of around 100 beats per minute at rest). For a totally stochastic heartbeat time series with averaged heart rate of 100 beats per minute, the average period of IMF 1 is 1.686 s. Then, we found that the intrinsic timescale of IMF 1 for a healthy sub-
494
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
Fig. 4. Analysis results of new scaling factor of EMD-based analysis and DFA scaling exponent. Pearson’s correlation coefficient between results using two different parameters is 0.799: (a) new scaling factor of EMD-based analysis; (b) DFA scaling exponent.
Fig. 5. Human heartbeat time series in breath-control experiment.
ject has a rhythm similar to respiration (between 3 and 5 s) but the timescale for a subject with CHF or AF has rhythm similar to IMF 1 decomposed from a totally stochastic time series (close to 1.68 s). Finally, according to the results shown in Table 1, the statistical difference between intrinsic timescales of healthy group (i.e., young or elderly) and group with heart disease (i.e., CHF or AF) is significant (P-value < 10−4 by t-test). But, the statistical difference between groups young and elderly is insignificant (P-value > 0.05 by t-test). Also, the statistical difference between groups CHF and AF is insignificant (P-value > 0.05 by t-test). Therefore, an intrinsic timescale of IMF 1 decomposed from a human heartbeat time series can be a clue to evaluate the RSA modulation of a cardiac system. In a healthy subject, the RSA modulation is clear and reflected by the intrinsic timescale of IMF 1. But, the RSA modulation is diminished and approaches to the characteristic
component decomposed from a stochastic time series for a subject with CHF or AF. Based on this finding, an intrinsic timescale of IMF 1 can be used as an indicator for the orderliness of a human heartbeat time series. A low value of intrinsic timescale similar to IMF 1 decomposed from a stochastic time series reflects a high degree of disorderliness, as characteristic of a subject with CHF or AF, whereas a high value indicates RSA modulation of a healthy subject. In addition, the receiver operating characteristic (ROC) curve can be used to examine the sensitivity and specificity of classification using different parameters. Therefore, the ROC curve was used to estimate the sensitivity and specificity of classification of subjects with or without heart disease. This used the both fractal property and intrinsic timescale of IMF 1 to show the critical parameter for classification. Fig. 6 shows that the intrinsic timescale of IMF 1 is a better indicator for distinguishing healthy subjects from subjects
Table 1 Statistical values for intrinsic timescale of IMF 1 for heartbeat intervals of four groups (i.e., healthy young, healthy elderly, CHF and AF). Group
Number of cases
Age
Young Elderly CHF AF
20 20 43 9
25.95 74.55 55.37 43.48
± ± ± ±
4.31 4.45 11.32 12.25
Severity
Intrinsic timescale of IMF 1 (s)
N/A N/A NYHA classes I, II, & III N/A
3.34 3.28 1.98 1.77
± ± ± ±
0.45 0.55 0.31 0.12
Notation Similar to the rhythm of respiration Similar to the rhythm of respiration Similar to IMF 1 decomposed from a random time series Similar to IMF 1 decomposed from a random time series
N/A means not available and NYHA is the abbreviation of New York Heart Association. Here, NYHA means the classification of New York Heart Association.
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
495
inverse relationship with aging showing that averaged amplitude of IMF 1 has a negative connection with aging. However, since the aging process is exceedingly complex, many variables should be taken into account in the experimental protocol for a further study, which aims to understand the underlying mechanisms of aging. Fig. 7 presents the statistical results of intrinsic timescale and averaged amplitude for four groups. Here, intrinsic timescales of healthy groups (i.e., young and elderly) are significantly different from those of groups with heart disease and averaged amplitudes of groups of healthy young and AF are significantly different from those of groups of elderly and CHF. However, IMF 1 reflects different underlying modulations for different groups. For healthy groups, it reflects the strength of RSA modulation. For example, the cardiac system of a healthy young subject has a stronger RSA modulation than that of a healthy elderly subject. The averaged amplitude of IMF 1 for a young subject is larger than that for an elderly subject. Therefore, the average amplitude of IMF 1 has a positive correlation with the strength of RSA modulation and a negative correlation with aging. Fig. 6. ROC curves for classification of subjects with or without heart disease using two different parameters of the fractal property (i.e., new scaling factor of EMDbased analysis) and averaged period of IMF 1.
with CHF or AF. As the physiological meaning, RSA modulation is an important characteristic of human cardiac system. Thus, intrinsic timescale of IMF 1 is an effective indicator in diagnosis of heart disease (i.e., CHF and AF in this study). 5.2. Connection between aging and RSA modulation As mentioned above, IMF 1 of human heartbeat time series for a healthy subject reflects RSA modulation of the cardiac system. An intrinsic timescale of IMF 1 was used to check with the rhythm of respiration as an indicator for subjects with or without heart disease. Furthermore, to understand the connection between aging and RSA modulation, we conducted a comparison between averaged amplitudes of IMF 1 for groups of healthy young and healthy elderly. These two groups have a significantly statistical difference in average amplitude of IMF 1 (P-value = 0.0039 by t-test), with the average amplitude of IMF 1 for the young group significantly bigger than that for the elderly group. Thus, RSA modulation has an
6. Discussion and conclusions In this study, we present a new approach for human heartbeat analysis, using an EMD-based analysis algorithm. We extended a finding based on the function of EMD, in which the hierarchy of the extracted mode functions may be exploited for obtaining the fractal property of time series, to develop a new scaling parameter to quantify the long-term fractal property of human heartbeat time series. According to the results of a numerical experiment, the new parameter had been proven to perform as a more accurate scaling index than the Hurst exponent by R/S analysis via the simulated time series of fractal Gaussian noise with various Hurst exponents from 0.1 to 0.9. Moreover, its results are highly consistent with those using the long-term scaling exponent by DFA in human heartbeat analysis. Theoretically, fractal property reflects the characteristic of selfsimilarity for a signal among continuous observing scales. For a stochastic time series, both DFA and the EMD-based analysis reflect a consistent characteristic of self-similarity by straight lines of logarithmic energy/timescale plots. However, a human heartbeat time series is not a totally stochastic time series, which contains dominant components caused by underling modulations. Hence, the
Fig. 7. Statistical values of intrinsic timescales and averaged amplitudes of IMF 1 for four groups.
496
J.-R. Yeh et al. / Medical Engineering & Physics 32 (2010) 490–496
logarithmic period-energy plots of EMD-based analysis are not a well-defined straight line for heartbeat time series of healthy subjects and some CHF patients. On contrast to EMD-based analysis, the trends are approximated by least square lines in DFA, which are insignificantly affected by internal modulation. Therefore, the analysis results of EMD-based analysis are more variable than those of DFA. Excluding concerns of the internal modulation in human heartbeat time series, the fractal property of human heartbeat time series is an important aspect for system identification of human cardiac system, consistent with the results of previous studies of human heartbeat analysis. Fractal property performs well to show the characteristic property of a human’s heartbeat time series for verifying an ill subject with CHF or AF, but it does not provide a critical parameter in system verification between subjects with or without heart disease. In addition, the statistical difference between the fractal properties of groups of young and elderly is insignificants. On the other hand, previous studies also showed that shortterm scaling exponent (˛1 ) of DFA works to quantify the ensemble responses of RSA and ANS modulations [21,22]. However, the application of short-term scaling exponent of DFA in human heartbeat analysis has some limitations. A basic limitation is that ˛1 works only for healthy subjects (i.e., subjects without heart disease). Thus, EMD had been used to investigate the intrinsic modulation in human heartbeat time series. According to our findings in this EMDbased analysis, IMF 1 reflects the RSA modulation of cardiac systems only for healthy subjects but not for subjects with CHF or AF. For subjects with CHF or AF, IMF 1 has similarly characteristic timescale to that of IMF 1 decomposed from a stochastic time series. This shows that the human heartbeat time series of subjects with CHF or AF contains high-frequency components, which diminishes the component of RSA modulation. This causes that short-term scaling exponent of DFA works well only for quantifying the RSA and ANS modulation of healthy subjects but not for subjects with CHF or AF. Furthermore, according to our new findings, the intrinsic timescale of IMF 1 can be used as an index of orderliness for human heartbeat time series. A value of intrinsic timescale of IMF 1 similar to the rhythm of respiration indicates the underlying orderliness of RSA modulation to the human cardiac system. On contrast to the orderliness caused by RSA modulation, IMF 1 with an intrinsic timescale of 1.68 s presents a high degree of randomness for human heartbeat time series. This finding is the most important contribution of our study. Simply speaking, RSA modulation is a key factor for diagnosis of subjects with or without heart disease (i.e., CHF or AF in this study) and the intrinsic timescale of IMF 1 decomposed from human heartbeat time series is a key indicator for RSA modulation. Moreover, average amplitude of IMF 1 for subjects with heart disease reflects the effects caused by the presence of noise. The present noise in human heartbeat time could be met in the high-frequency part and embedded in IMF 1. The noise present in heartbeat time series is comparatively weak to result small amplitude of IMF 1 for a CHF patient. But, the present noise dominates the heartbeat time series and produces large amplitude of IMF 1 for an AF patient. Therefore, amplitude of IMF 1 reflects the strength of noise presenting in heartbeat time series as a key parameter for distinguishing the subject with CHF or AF. In summary, our previous study showed there is another dimension, different from the fractal property evaluated by the original scaling exponent of DFA, which is a new parameter for the distribution of intrinsic timescales. In this study, we found that RSA contributes to the fundamental modulation of the human cardiac system in healthy subjects. Thus, irregularity of RSA modulation,
as shown in IMF 1, is a significant sign to a cardiac system with CHF or AF. Furthermore, EMD is a powerful tool to extract the RSA modulation from human heartbeat time series, which is the most important contribution in the EMD-based analysis. Acknowledgements The authors wish to thank Prof. Ary L. Goldberger and Dr. C.K. Peng (the director and co-director of the Rey Institute for Nonlinear Dynamics in Medicine at the Beth Israel Deaconess Medical Center of Harvard Medical School) for valuable discussion. We gratefully acknowledge support from the National Science Council (NSC) of Taiwan (Grant Number NSC96-2221-E-155-015-MY3-2) for supporting this research. Conflict of interest statement The authors want to disclose there are no potential conflicts of interest with other people or organizations in this study. References [1] Costa M, Goldberger AL, Peng CK. Multiscale entropy analysis of complex physiologic time series. Phys Rev Lett 2002;89(6):068102. [2] Goldberger AL, Peng CK, Lipsitz LA. What is physiologic complexity and how does it change with aging and disease? Neurobiol Aging 2002;23:23–6. [3] Ho KK, Moody GB, Peng CK, Mietus JE, Larson MG, Levy D, et al. Predicting survival in heart failure case and control subjects by use of fully automated methods for deriving nonlinear and conventional indices of heart rate dynamics. Circulation 1997;96:842–8. [4] Francis DP, Willson K, Georgiadou P, Wensel R, Davies LC, Coats A, et al. Quantitative general theory for periodic breathing in chronic heart failure and its clinical implications. Circulation 2000;102:2214–21. [5] Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc Roy Soc Lond A 1998;454:903–95. [6] Wu Z, Huang NE, Long ST, Peng CK. On the trend, detrending, and variability of nonlinear and nonstationary time series. PNAS 2007;104(38):14889–94. [7] Hadjileontiadis LJ. A novel technique for denoising explosive lung sounds empirical mode decomposition and fractal dimension filter. IEEE Eng Med Biol Mag 2007;26(1):30–9. [8] Wu Z, Huang NE. A study of the characteristics of white noise using the empirical mode decomposition method. Proc Roy Soc Lond A 2004;460:1597–611. [9] Flandrin P, Rilling G, Goncalces P. Empirical mode decomposition as a filter bank. IEEE Signal Process Lett 2004;11:112–4. [10] Yeh JR, Fan SZ, Shieh JS. Human heart beat analysis using a modified algorithm of detrended fluctuation analysis based on empirical mode decomposition. Med Eng Phys 2009;31:92–100. [11] Hirsch JA, Bishop B. Respiratory sinus arrhythmia in humans: how breathing pattern modulates heart rate. Am J Physiol Heart Circ Physio 1981;241:H620–9. [12] Berntson GG, Carcioppo JT, Quigley KS. Cardiac psychophysiology and autonomic pace in humans: empirical perspectives and conceptual implications. Psychol Bull 1993;114(2):296–322. [13] Craft N, Schwartz JB. Effects of age on intrinsic heart rate, heart rate variability, and AV conduction in healthy humans. Am J Physiol 1995;268:H1441–52. [14] Databases are available at, http://www.physionet.org/, See Goldberger AL, et al. Circulation 2000;101:E215. [15] Flandrin P, Goncalves P. Empirical mode decomposition as a data-driven wavelet-like expansions. Int J Wavelet Multires Info Proc 2004;2:477–96. [16] Mandelbrot BB, van Ness JW. Fractional Brownian motions, fractional noises and applications. SIAM Rev 1968;10:422–37. [17] Hurst HE. Long-term storage capacity of reservoirs. Trans Am Soc Civil Eng 1951;116:770. [18] Wood AT, Chan G. Simulation of stationary processes in [0,1]d . J Comput Graph Stat 1994;3:409–32. [19] Alessio E, Carbone A, Castelli G, Frappiettro V. Second-order moving average and scaling of stochastic time series. Eur Phys J B 2002;27:197–200. [20] Glantz SA. Primer of biostatistics. 6th ed. Singapore: McGraw-Hill Inc.; 2005. [21] Peng CK, Havlin S, Stanley HE, Goldberger AL. Quantification of scaling exponents and crossover phenomena in nonstationary heart-beat time series. Chaos 1995;5:82–7. [22] Mourot L, Boudaddi M, Gandelin E, Cappelle S, Nguyen NU, Wolf JP, et al. Conditions of autonomic reciprocal interplay versus autonomic co-activation: effects on non-linear heart rate dynamics. Autonom Neurosci: Basic Clin 2007;137:27–36.