Chaos-Based Analysis of Heart Rate Variability Time Series in Obstructive Sleep Apnea Subjects.

Obstructive sleep apnea (OSA) is a common disorder which can cause periodic fluctuations in heart rate. To diagnose sleep apnea, some studies analyze electrocardiogram (ECG) signals by adopting chaos-based analysis. This research is going to specifically focus on whether it is possible to use chaos-based analysis of heart rate variability (HRV) signals rather than using chaotic analysis of ECG signals to diagnose OSA. While conventional studies mostly use chaos-based analysis of ECG signals to detect OSA, here, we apply correlation dimension (CD) as a chaotic index to analyze HRV data in OSA patients. For this purpose, 17 patients with OSA and 9 healthy individuals referred to a sleep clinic in Isfahan/Iran are studied, and their HRV time series were extracted from 1-h ECG signals recorded overnight. The preliminary step to calculate CD is phase-space reconstruction of the system based on HRV time series. Corresponding parameters, including embedding dimension and lag time, are estimated optimally using enhanced related methods, and then CD is calculated using Grassberger-Procaccia algorithm. Moreover, to evaluate our results, detrended fluctuation analysis (DFA), one of the well-known nonlinear methods in HRV analysis to detect OSA, is also applied to our data and the result is compared with those obtained from CD analysis of HRV. CD index with P < 0.005 indicates a significant difference in nonlinear dynamics of HRV signals detected from OSA patients and healthy individuals. Copyright:

Introduction

Sleep apnea or pauses in breathing for more than 10 s followed by awakening, is one of the common sleep disorders. Sleep apnea is divided into three categories: obstructive sleep apnea (OSA), central sleep apnea, and mixed (complex) sleep apnea. OSA, the most common one, is caused by a local blockage or obstruction of the upper airway and is usually accompanied by exhale and deep inhale at the end of blockage. Obesity, aging, narrowed airways, smoking, and alcohol consumption are considered as the certain risk factors of OSA, and men are more likely to have sleep apnea.[1] Sleep apnea reduces the quality of life since it may cause depression, irritation, and lack of sleep. Investigations reveal that OSA is the reason of many car accidents and work events;[23] furthermore, according to the prevalence of cardiovascular diseases among those suffering from sleep apnea, many studies have been done to explore the relation of hypertension and heart diseases to sleep apnea. Studies have shown that untreated OSA can lead to high blood pressure, cardiac arrhythmias, heart failure, and brain stroke as well.[456] Studies on OSA reveal that apnea is accompanied by periodic variations in heart rate.[7] Irregular heartbeats appear in the form of bradycardia (decrease in heart rate to <60 times/min) during apnea, and in the form of tachycardia (increase in heart rate to more than 100 times/min) after that; moreover, according to the research done on sympathetic neural activity, sympathetic activation increases during apnea.[8] In patients with OSA, sympathetic activation continues during the day that can increase the risk of getting cardiovascular diseases.[9] There have always been many interests in finding simple and accurate methods for the study of or detection of OSA; accordingly, some research on detection of sleep apnea has been already done using different biological signals such as electroencephalography signal[10] and electrocardiogram (ECG) signal with or without respiratory signal.[1112] Furthermore, blood oxygen measurements[13] and feature assessments of acoustic signal produced by snoring[14] have been already used for sleep apnea detection. According to the studies,[7] heart rate variability (HRV), variation in the time intervals between successive heartbeats, has been also used efficiently to detect sleep apnea. HRV has become a conventional marker to detect different failures in autonomous nervous and cardiovascular system.[15] HRV signals are typically extracted from ECG signals by measuring the R–R intervals. Like other physiological signals, HRV can be analyzed using time/frequency-domain, nonlinear, and chaotic methods. Time-domain methods calculate the statistical features of R–R intervals. They can only provide the simplest insight into the cardiovascular system. Frequency-domain indexes, which are mostly based on spectral bands, can reveal the autonomic nervous system activities. However, nonlinear methods can describe hidden dynamics and some nonstationary conditions of the cardiovascular system.[16] The human heart acts like a nonlinear oscillator, and HRV is one observed variable from a complex system (cardiovascular system); therefore, complexities in the cardiovascular system can be revealed using nonlinear analysis of HRV.[1718] Based on the study done by Voss et al. [19] due to nonlinear nature of biological signals, nonlinear techniques can be efficiently used for the analysis of these signals; moreover, in biological signal analysis, nonlinear methods can provide complementary information along with linear methods (time/frequency-domain methods). The conventional nonlinear methods used to analyze HRV in sleep studies include detrended fluctuation analysis (DFA), approximate entropy (ApEn), and sample entropy (SampEn) as well as chaotic methods such as recurrence quantification analysis (RQA), correlation dimension (CD), and Lyapunov exponent (LE). Some sleep studies focused on comparing the use of nonlinear methods with linear techniques for analysis of HRV in sleep apnea patients. In a study performed by Penzel et al.,[20] frequency-domain analysis and DFA were used to investigate the effect of sleep apnea on HRV. Their results suggest that in comparison with spectral analysis, DFA as a nonlinear method can better depict fluctuations caused by sleep apnea in HRV signal. In another study, Al-Angari and Sahakian[21] have applied SampEn as a nonlinear method to assess the complexity of HRV behavior due to OSA. According to their results, in healthy individuals, HRV pattern is considerably more complex than in patients with OSA. They also compared the results of the spectral analysis as a linear method with those obtained from SampEn. They claimed that although both methods are accurate, SampEn calculations are simpler in comparison with spectral methods.

Some noticeable nonlinear methods in sleep research are chaos-based analysis. CD is a well-known chaotic index used in some sleep apnea studies. In a study, Miyata et al.[22] used CD to investigate nonlinear features and chaotic behavior of respiratory movement and breath to breath fluctuations in patients with OSA/hypopnea syndrome during wakefulness with eyes closed. They claimed that applying CD to respiratory movement can be used to efficiently identify patients with this syndrome during wakefulness. In their previous study,[23] they also showed using chaotic analysis of respiratory movement cannot be useful to detect OSA/hypopnea syndrome during apneic sleep. In a recent work, Moeynoi and Kitjaidure[24] explored RQA, CD, and DFA to extract nonlinear features from HRV signals and ECG derived respiratory signals to detect sleep apnea. However, their main focus of this work is on introducing canonical correlation analysis as a method for dimensional reduction of sleep apnea features extracted from ECG for better classification. In such works, discrimination abilities of chaotic methods have been rarely investigated in a specific way. These methods have been only used for classification purposes. In one study focused on chaotic analysis, Acharya et al.[25] applied nonlinear parameters including ApEn, fractal dimension, CD, largest LE, and Hurst exponent to ECG signals to distinguish sleep apnea patients. Nevertheless, in this study, phase-space reconstruction which is the primary step to calculate CD was performed using ECG signals as time series. In another study, Jafari[26] used nonlinear methods consisting of DFA, CD, large LE, and spectral entropy to detect OSA. He found that detrended fluctuation and CD are more significant for OSA detection among other methods he applied, but like the work done by Acharya et al.,[25] he used ECG signals for phase-space reconstruction. LE is another well-known chaotic index with the ability to detect chaotic behavior of biological signals. Jafari showed large LE is also a significant index for OSA detection. However, Zapanta et al.[27] introduced a new method called numerical titration technique instead of LE to detect chaos in HRV of children with OSA. They claimed that numerical titration technique could provide robust support for the presence of chaotic behavior in HRV. In spite of such studies worked on chaotic analysis in sleep research, up to our knowledge, sleep apnea studies have rarely focused on calculation chaotic indexes using HRV time series. As above mentioned, in the study by Jafari,[26] ECG signals are used as time series for the phase-space reconstruction as well as estimation of embedding dimension and lag time parameters, which are necessary parameters to calculate chaotic indexes; hence, this question arises that whether it is possible to use HRV of OSA patients rather than their ECG signals as time series to calculate chaotic indexes or not. Thus, the main aim of this study is specific concentration on the investigation of CD calculation using HRV as time series. In fact, while conventional studies such as the study by Acharya et al.[25] or by Jafari,[26] mostly use chaos-based analysis of ECG signals to study OSA, in this research, we concentrate to calculate CD using HRV signals instead of using ECG signals. We are going to find out whether there is a significant difference in nonlinear dynamics of HRV signals between OSA patients and healthy individuals. We are also going to compare our results of CD calculated by HRV data as time series with the results obtained from DFA, which is one of the most commonly used nonlinear methods in sleep apnea studies.[20]

Materials and Methods

Data

In this research, the exploited data set is obtained from a cross-sectional study[28] on patients suffering from OSA. This data set contains both patients and normal individuals. The patient group consisted of individuals with untreated OSA who were referred to the Bamdad Sleep Clinic in Isfahan, Iran, and were diagnosed as mild-to-severe OSA, based on Apnea-Hypopnea Index (≥5 events/h), and the age- and sex-matched control group consisted of healthy individuals without any sleep or cardiac complaints. All patients and normal individuals were visited and examined by a cardiologist, and those with a history or clinical evidence of heart failure, ischemic heart disease, cardiomyopathy, valvular heart disease, arrhythmia, persistent atrial fibrillation, bundle branch block, pericarditis, electrolyte abnormalities, renal failure, pulmonary disease, thyroid dysfunction, hypertension, diabetes, and use of medications affecting the electrocardiographic parameters, were not included into the study. The study protocol was approved by the Ethics Committee of the Isfahan University of Medical Sciences (Grant # 392101), and informed consent was provided by all patients and normal individuals. Patients had undergone an overnight polysomnographic study at the sleep clinic while all of them had effective sleep.

In this study, a 24-h electrocardiography was performed for patients and normal individuals by using a digital Holter recorder (H200 recorder, Kavoshgaran Teb Kharazmi Co, Iran), with six channels (three analogs I, II, and III and three digital, aVR, aVL, and aVF channels). Software was developed by the Biomedical Engineering Department of the University of Isfahan to measure ECG parameters. This software detects Q, R, S, P, and T points using gradient-based methods; then calculates important time intervals in ECG signal. RR intervals of ECGs achieved by H200 recorder for both patient and normal groups were used as an HRV data set in this research. The HRV time series extracted from 1-h ECG signals (recorded during the night) related to 17 patients with OSA and 9 healthy controls are analyzed.[7]

Preprocessing

To analyze HRV related to sleep apnea, Penzel et al.[20] have selected an experimental feature for automated preprocessing of HRV time series or RR intervals. The feature is as follows:

0.33 s

The interval between two successive RR intervals is just 0.66 s.

If >1% of RR time intervals of a signal do not meet the above conditions, the signal cannot be used for future analysis. Of course, those HRV samples that exceed the above limit have to be discarded from the accepted signals.

Since the resulted signals were too noisy, after the evaluation of the above conditions, a median filter was applied to the signal to preserve the latent dynamics in noisy data.

Correlation dimension

CD is a measure of the complexity of the system that shows variability and irregularity of the process. For calculation of CD as a chaotic index, reconstruction of a phase-space equivalent to the original phase space in topological aspects is the first important step. Takens’ theory is often used for phase-space reconstruction,[29] that is called as the method of delays. In this method, for a time series derived from a deterministic dynamical system, the delay coordinates are used to form an m-dimensional vector space as follows:

x(t)=f(x(t–τ),x(t–2τ),......,x(t–(m-1)τ))      (1)

Where scalar m called embedding dimension and scalar τ called time delay.

This theory only explains the possibility of phase-space reconstruction with two parameters, embedding dimension and time delay, but it does not provide any information about the way of acquiring these parameters. Time delay or τ expresses the distance between components of the delay vector and embedding dimension or m expresses the components in each delay vector.

The success of CD analysis depends on suitable reconstruction of phase space that is based on the calculation of embedding dimension and lag time. There are several ways to obtain these parameters. In this study, the mutual information method is selected to achieve the optimized value for time delay;[30] Indeed, the mutual information (M (τ)) is obtained for different τ values from (2).

Where P is the probability of finding a data point in the distance subscripted by i and P is the joint probability of an observation in the distance subscripted by i while the next observation lies in the distance subscripted by j after the time delay τ. The first minimum of the M (τ) is considered as the optimum time delay for phase-space reconstruction. Figure 1 shows M (τ) calculated for several τ values for the HRV signal of a patient with OSA.

Figure 1

Mutual information diagram for a heart rate variability signal from a patient with obstructive sleep apnea

For determination of the optimized embedding dimension, the false nearest neighbor method is exploited which is one of the most precise methods for the selection of optimum embedding dimension in reconstructing a low-dimension phase space.[31] In this method, optimum embedding dimension is selected as the first dimension for which the amount on nearest neighbors equals or approaches to zero.

After reconstruction of phase space, CD can be calculated using the Grassberger–Procaccia algorithm which in this regard is one of the well-known and optimized approaches. Relation or correlation between delay vectors extracted from phase-space reconstruction is measured with correlation integral. Correlation integral is an estimation of the probability of two delay vectors being at a distance less than r. Correlation integral is resulted from (3).

Where Θ shows Heaviside function, N is the number of points of the time series, and r is the distance under study; x and x refer to the points located on the trajectories of the reconstructed phase space. Heaviside function is shown in (4).

Then, CD is calculated using (5).

When log C(r) is plotted as a function of log r, three distinct regions appear: Null region, the region in which we could not find two points, located at a distance less than r; Saturation region, the region in which C(r) does not grow by increasing r; and Scaling region, the region between null and saturation regions. CD equals to the slope of this curve in the third region.[3233] Figure 2 shows log C (r) as a function of log r extracted from the HRV signal related to a patient with OSA.

Figure 2

The diagram of log C (r) as a function of log r for a heart rate variability signal from a patient with obstructive sleep

Detrended fluctuation analysis

DFA is used to determine the self-similarity and fractal structure in a time-series. In DFA, root mean square (RMS) fluctuation resulted from a detrended and integrated time series is measured using different frames. These values are plotted as a function of applied frames in a log-log plot. First, the time series y(k) with the total length of N, is integrated with the length of K ∈ N depicted in (6).

Where y(i) is data point number i from y(k) time series and y is the mean value of the time series. The integration process reveals that the main signal is nonstationary. After integration, the integrated time series is equally divided into subsets; each subset with length L. Then, a line obtained using the least square error method fits to data of each section and shows the trend of that section. After that, local trends related to each section are removed from the integrated time series Y (k). For higher-order trend removal, a higher-order polynomial should fit to the time series. At last, RMS of the final time series is calculated as in (7) and (8).

Where Yv(k) shows lines fitted to each section (with length L) of the time series Y(k), N shows the number of sections after division of time series into sections with length L. The above calculations should be continued until a meaningful relevance between F(L) and L is defined. F(L) typically grows in proportion with L. Slope of log F(L) as a function of log L achieves the self-similarity parameter called α, which is the measurable parameter in DFA.[34] α is calculated using (9).

The primary step toward calculations of CD is phase-space reconstruction of the system based on HRV time series. Two corresponding important parameters, including embedding dimension and lag time, are estimated optimally using enhanced related methods. The time delay is acquired through mutual information method. The first minimum of the mutual information diagram is considered as the optimum time delay for the phase-space reconstruction. M (τ) is calculated for several τ values for the HRV signal of a patient with OSA and is shown in Figure 1. Time delay for all the HRV signals are calculated and the dominant value (τ = 20) is selected for the phase-space reconstruction. Then, embedding dimension is achieved using the nearest neighbor method. The first dimension, for which the number of nearest neighbors approaches to 0, is selected as the embedding dimension. After the evaluation of embedding dimensions obtained from all HRV signals, the dominant value (m = 5) is used for the phase-space reconstruction; therefore, CD is calculated using τ = 20 as the time delay and m = 5 as the embedding dimension. Table 1 shows CDs calculated for HRV signals recorded from both patients with OSA and healthy individuals. The results of t-test analysis performed on CD values are statistically significant.[35] The CD with P < 0.005 indicates that there is significant difference in nonlinear dynamics of HRV signals related to OSA patients and healthy individuals; indeed, our result indicates the possibility of using HRV data in sleep apnea as time series to calculate chaotic indexes.

Table 1

Results

Features	Mean±SD		P
	Normal	OSA
CD	2.696±0.194	2.933±0.184	0.003
α	1.159±0.11	1.063±0.132	0.044

OSA – Obstructive sleep apnea; CD – Correlation dimension

To evaluate our results with another method, we also apply DFA to our data. In DFA, RMS should be repeatedly calculated for at least 30 frames. In this study, according to the work done by Penzel et al.,[20] 60 frames between L = 4 and L = 1000 are taken into account. Table 1 shows α values obtained from DFA on HRV signals for patients with OSA and healthy controls. As mentioned above, α is the measurable parameter resulted from DFA. The results of t-test analysis performed on α values are statistically significant. Our obtained results shown in Table 1 indicate that in comparison with DFA, CD as a chaotic method which is calculated using HRV time series can better quantify HRV changes.

Discussion and Conclusion

Chaos-based methods can assist to realize hidden dynamical features of biological time series, so they can be used to detect serious disease such as sleep apnea. In most chaotic analysis methods such as CD, phase-space reconstruction is the first step. Currently, most of the sleep apnea studies calculating CD exploit ECG signals as time series for phase-space reconstruction. Up to our knowledge, studies on sleep apnea disorder rarely focused on independent investigation of the use of HRV as time series to perform phase-space reconstruction and to calculate CD. To this end, we attempted to use HRV data as time series to calculate CD. Statistical tests such as t-test typically are necessary to evaluate the discrimination ability of every feature extracted from time series analysis. Our results of using t-test on CD obtained from HRV time series with P = 0.003 show a significant difference in nonlinear dynamics of HRV signals detected from OSA patients and healthy individuals, so this method has a very significant ability to distinguish OSA patients from normal people. A few studies on OSA using chaotic analysis can obtain practically significant results; for instance, Miyata et al.[22] obtained results with P < 0.01 by applying CD to respiratory movement for diagnosis of OSA/hypopnea syndrome during wakefulness, or Zapanta et al.[27] got results with P < 0.05 to detect chaos in HRV of children with OSA. Faust et al.[36] in their study claim that from research experience, features with low P value are more practical to discriminate patients from healthy people. Based on their claim, having features with low P value is necessary to design an automatic diagnosis system. Therefore, our statistical results with practical significance prove that the CD obtained from HRV time series can be used to explain an appropriate classifier for automatic detection of OSA. Design of a classifier based on distinguish ability of our method is another work that we already aim to propose. Moreover, we also used DFA (a commonly used method for OSA detection) to examine our HRV data. Comparing our obtained results from CD method with our results from DFA also indicates that HRV changes of sleep apnea patients can be efficiently quantified using CD. As shown in Table 1, our findings indicate more discrimination ability for CD than DFA. In conclusion, we attempted to specifically investigate HRV as time series to perform the phase-space reconstruction and to calculate CD as well, and our results as discussed above can be efficiently useful in study of OSA.

Financial support and sponsorship

None.

Conflicts of interest

There are no conflicts of interest.

BIOGRAPHIES

Shiva Naghsh was born in 1989. She received her B.Sc. degree in biomedical engineering in 2011 and the M.Sc. degree in electrical engineering in 2014 from the University of Isfahan, Isfahan, Iran. She is currently working as an R&D engineer. Her research interests include biological signal processing, time series analysis, and data analysis.

Email: sh.naghsh@gmail.com

Mohammad Ataei received the B.Sc. degree from the Isfahan University of Technology, Iran, in 1994, the M.Sc. degree from the Iran University of Science and Technology, Iran, in 1997, and the Ph.D. degree from the K. N. Toosi University of Technology, Iran, in 2004 (joint project with the University of Bremen, Germany- DAAD scholarship), all in control Engineering. He is currently a full Professor with the Department of Electrical Engineering, University of Isfahan, Iran. His main areas of research interest are control theory and applications, nonlinear control, and chaotic systems’ analysis and control.

Email: ataei@eng.ui.ac.ir

Mohammadreza Yazdchi obtained his B.Sc. in Electrical Engineering from the Isfahan University of Technology in 1997. He received his M.Sc. And Ph.D. in Biomedical Engineering from Amirkabir University of Technology in 2000 and 2006 respectively. He is currently an Associate Professor in the Department of Biomedical Engineering at the University of Isfahan since 2007. His research interests include Biomedical Signal Processing, Medical Image Processing ,Biomedical Instrumentation, Biologically Inspired Computing and Bio-Inspired Engineering.

Email: yazdchi@eng.ui.ac.ir

Mohammad Hashemi was born in Isfahan in 1965. He studied medicine at Isfahan University of Medical Sciences (IUMS) and graduated in 1990. He continued his studies to become specialized in Cardiology in IUMS until 1994 and begin fellowship of Interventional Cardiology in 1997 in Iran University Of Medical Sciences. In 1999, he joined the department of Internal Medicine of IUMS until now. He is currently a full professor of Cardiology. He is interested on Endothelial Function as academic research.

Email: mohammad_hashemi@med.mui.ac.ir