Structures of the recurrence plot of heart rate variability signal as a tool for predicting the onset of paroxysmal atrial fibrillation.

This paper aims to propose an effective paroxysmal atrial fibrillation (PAF) predictor which is based on the analysis of the heart rate variability (HRV) signal. Predicting the onset of PAF, based on non-invasive techniques, is clinically important and can be invaluable in order to avoid useless therapeutic interventions and to minimize the risks for the patients. This method consists of four steps: Preprocessing, feature extraction, feature reduction, and classification. In the first step, the QRS complexes are detected from the electrocardiogram (ECG) signal and then the HRV signal is extracted. In the next step, the recurrence plot (RP) of HRV signal is obtained and six features are extracted to characterize the basic patterns of the RP. These features consist of length of longest diagonal segments, average length of the diagonal lines, entropy, trapping time, length of longest vertical line, and recurrence trend. In the third step, these features are reduced to three features by the linear discriminant analysis (LDA) technique. Using LDA not only reduces the number of the input features, but also increases the classification accuracy by selecting the most discriminating features. Finally, a support vector machine-based classifier is used to classify the HRV signals. The performance of the proposed method in prediction of PAF episodes was evaluated using the Atrial Fibrillation Prediction Database which consists of both 30-minutes ECG recordings end just prior to the onset of PAF and segments at least 45 min distant from any PAF events. The obtained sensitivity, specificity, and positive predictivity were 96.55%, 100%, and 100%, respectively.

INTRODUCTION

Atrial fibrillation (AF) is the most common cardiac arrhythmia and entails an increased risk of thromboembolic events.[1] Although AF is not considered a life-threatening arrhythmia, it may severely impact the quality of life and increase the risk of stroke. About 15% of strokes occur in people with AF.[2]

Prevention and treatment of AF is still far from satisfactory. The aim of the therapy is to prevent stroke and regain sinus rhythm.[3] Clinically, AF presents itself in different forms. Often it starts as paroxysmal (self-terminating) and becomes more persistent with time. Paroxysmal AF (PAF) is defined as attacks of AF lasting from two minutes to less than seven days and spontaneously revert to normal sinus rhythm. Permanent AF (non-terminating) is defined as lasting more than seven days and sinus rhythm cannot be restored or maintained. About 18% of PAF evolve to permanent AF over four years.[3]

Development of accurate predictors of the onset of PAF is clinically important because of the increasing possibility of electrically stabilizing and preventing the onset of atrial arrhythmias with different atrial pacing techniques.[4] The maintenance of sinus rhythm can lead to decreased symptoms and possibly a decrease in the atrial remodeling that causes increased susceptibility to future episodes of PAF.[5] In addition, there may be a reduction in the risk of strokes and thromboembolic events.

In the last decades, several methods have focused on finding algorithms able to predict PAF by the analysis of surface electrocardiographic (ECG) records. Most of these methods have been evaluated on 30-min ECG segments from Atrial Fibrillation Prediction Database (AFPDB) and the concept of PAF prediction in the literature, as our proposed method, is to discriminate between ECG segments far from PAF and segments before PAF event. Zong et al.[6] studied the number and timing of premature atrial complexes (PACs) in each segment and found that not only the number of PACs increases in episodes preceding PAF, but also these complexes most occur near the end of these episodes. Thong et al.[7] developed an algorithm based on the number of isolated PACs not followed by a regular RR-interval, runs of atrial bigeminy and trigeminy, and the length of any short run of paroxysmal atrial tachycardia. Their study showed that an increase in activity detected by any of these three criteria is an indication of an imminent episode of PAF. Vikman et al.[8] calculated the approximate entropy (ApEn) and short term scaling exponent, α1, of heart rate variability (HRV) over 20-min periods and found that a reduced complexity of HRV dynamics and altered fractal properties usually precede the onset of PAF as indicated by decreasing value for ApEn and α1. Chesnokov[9] combined complexity and spectral analysis of the 30-min HRV segment from AFPDB and found that there is statistically significant increase in the very low frequency (VLF) band, low frequency (LF) band, and high frequency (HF) band for the records immediately before PAF compared to distant ones, but LF/HF ratio doesn’t discriminate these two groups with statistical significance. He also showed that complexity features like ApEn, sample entropy (SmEn), and their multiscale versions exhibit smaller values in episodes preceding PAF compared to distant ones. In our previous work, we have proposed an algorithm for predicting the PAF attacks which was based on the spectrum, bispectrum, and nonlinear features of the HRV signal.[10] We showed that the frequency content and phase relations between frequency components of the HRV signal change before PAF events. Moreover, the nonlinear analysis can map the heart rate irregularities in the feature space and it leads to better understanding of the system dynamics before PAF attacks.

In recent years, a number of methods have been devised to compute dynamical features from time series.[11] Such features are the information dimension, entropy, Lyapunov exponents, dimension spectrum, etc. Recurrence is a fundamental property of dynamical systems, which can be exploited to characterize the system's behavior in phase space. A powerful tool for their visualization and analysis called recurrence plot (RP) was introduced in the late 1980's.[12] RP is the graphical representation of a binary symmetric square matrix which encodes the times when two states are in close proximity (i.e. neighbors in phase space). Based on such a recurrence matrix, a large and diverse amount of information on the dynamics of the system can be extracted and statistically quantified using recurrence quantification analysis (RQA).[13] No mathematical assumptions regarding the data and the generating systems constrain the construction of the RP, thus, this tool is particularly suitable for the analysis of physiological signals which are often non-stationary.[14] RPs have been applied in various biological system studies such as neuronal spike trains,[15] protein sequence,[16] analysis of EEG recordings,[1718] and classification of ECG abnormalities.[19]

In this paper, an accurate PAF prediction method is proposed, in which RQA is used to quantify the RP of the HRV signal. These features are based on the recurrence point density and the diagonal and vertical line structures of the RP. We have used six features to form an effective feature vector. Afterward, we have reduced these features to three new features by means of the LDA algorithm and fed them to a support vector machine (SVM)-based classifier to discriminate the episodes preceding a PAF event from the episodes being distant from any PAF event. SVM is a machine-learning technique which has established itself as a powerful tool in many classification problems. Simply stated, the SVM identifies the best separating hyperplane (the plane with maximum margins) between the two classes of the training samples within the feature space by focusing on the training cases placed at the edge of the class descriptors. In this way, not only an optimal hyperplane is fitted, but also less training samples are effectively used; thus, high classification accuracy is achieved with small training sets.[20] The results show that combining the extracted features and the SVM classifier yields an accurate predictor of PAF attacks.

In continue the details of the proposed algorithm for prediction of PAF is presented. Section 2 explains the database used for evaluation of the algorithm. Section 3 provides the overall block diagram of the proposed algorithm and the details of each block. Section 4 presents the results and discussion of the study. Finally, the conclusions are given in section 5.

DATABASE

The HRV data used in this work generated from the archived ECG signals provided by the Atrial Fibrillation Prediction Database (AFPDB).[21] This database consists of both 30-min ECG segments end just prior to the onset of PAF and segments at least 45 min distant from any PAF events. Each ECG record consists of two-channel traces from a Holter recording sampled at 128 Hz and 12-bit resolution. From this database, we have extracted 106 30-min ECG segments from 53 different patients. In other words, we have two ECG segments per patient: One segment contains the ECG immediately preceding an episode of PAF and the other segment contains the ECG during a period that is distant from any episode of PAF.

PROPOSED METHOD

Preprocessing

The aim of the preprocessing stage is to extract the HRV signals from the ECG signals within the database. In general, preprocessing can be affected by many interfering signals contaminating the ECG signal such as the 50 Hz power line interference, the interferences from EMG signals, and the baseline wandering. Therefore, in the preprocessing stage, these interfering noises are eliminated first by means of a 5-15 Hz band pass filter. Next, in order to detect the R peaks of the ECG, Hamilton and Tompkins algorithm[2223] is employed. The HRV signal is then constructed by measuring the time intervals between the successive R peaks.

Since ectopic beats and other artifacts are known to have serious impact on the results of the HRV analysis, we have edited the HRV signal to handle the ectopic beats and artifacts. To this end, we have used the method recommended by Kamath et al.[24] in which the HRV segments with duration of less than 80% of the duration of the previous normal beat have been regarded as ectopic beats. After detection of the ectopic beats, we have replaced the non-normal RR-intervals with new cubic spline interpolated RR-intervals.[25]

Feature Extraction

In this step, we attempt to extract the features from HRV signal which can be used as good markers for the prediction of PAF.

Recurrence Plot

The first step in the analysis of a signal using RP theory is the reconstruction of the phase space of the signal. A frequently used method for the reconstruction is the time delay method.[26] In this approach, the time series s (i=1, 2.,., N) of length N is embedded into a m-dimensional space with the time-delay (τ) technique

x = [s s ··· s] , i = 1,2,...,N – (m–1)τ.      (1)

For the analysis, both embedding parameters, the dimension m and the delay τ, have to be chosen appropriately. Different approaches for the estimation of the smallest sufficient embedding dimension (e.g. the false nearest-neighbors algorithm[27]), as well as for an appropriate time delay (e.g. the auto-correlation function, the mutual information function[28]) have been proposed. The false nearest-neighbors algorithm identifies the number of “false nearest neighbors” points that appear to be nearest neighbors because the embedding dimension is too small of every point in the phase space. When the number of false nearest neighbors drops to zero, we have embedded the time series to proper dimensional space.[27] For discontinuous signals such as RR-intervals extracted from the continuous ECG signals, the delay is best set to 1.[29] The explanations about setting the appropriate values for embedding parameters in our proposed method will be presented at the end of this section.

After the state space reconstruction, the RP of a signal can be obtained. In this step, the M×M recurrence matrix, the elements of which can be calculated using equation (2), is derived.

R = Θ(ε – ||x – x||),i,j = 1,2,…M,      (2)

where M=N-(m-1)τ, ε is a threshold distance,||·|| is the norm (e.g., the Euclidean norm) and Θ(x) is the Heaviside function. This means that if two phase space vectors x and x are sufficiently close together, then Ri,j = 1; otherwise Ri,j =0.

The RP is obtained by plotting the recurrence matrix. This means if the distance between x and x is less than ε, then a dot is placed at (i, j) in the RP. A crucial parameter of a RP is the threshold ε. Therefore, special attention has to be required for its choice. If ε is chosen too small, there may be almost no recurrence points and we cannot learn anything about the recurrence structure of the underlying system. On the other hand, if ε is chosen too large, almost every point is a neighbor of every other point, which leads to a lot of artifacts. A too large ε includes also points into the neighborhood which are simple consecutive points on the trajectory.[30] Several criteria for the choice of the cutoff distance ε have been proposed.[3132] One approach uses a fixed number of neighbors, Nn, for every point of the trajectory, called fixed amount of nearest neighbors (FAN).[12] In this approach, the cutoff distance, εi, changes for each state x to ensure that all columns of the RP have the same recurrence density. Using this neighborhood criterion, εi can be adjusted in such a way that the recurrence rate (REC) has a fixed predetermined value (i.e., REC=Nn/N).[30]

In this study, each 30-min HRV signal from AFPDB is split into five-min segments for the analysis, so there are six segments for each recording. Figures 1and 2 show the extracted RPs of the last 15-min of an episode far from PAF (record P05) and an episode before PAF (record P02), respectively. The optimum value for embedding dimension based on the criteria of the percentage of false nearest neighbors being less than 1%, is m=7 for the proper state space reconstruction of the HRV segments. The value of delay, τ, is set to one sample and no points in the HRV data are skipped for the state space reconstruction.

Figure 1

Recurrence plots of the last 15-min segment of the episode far from PAF (record P05). Forth to sixth 5-min segments of the episode are shown in (a) to (c)

Figure 2

Recurrence plots of the last 15-min segment of the episode before PAF (record P02). Fourth to sixth 5-min segments of the episode are shown in (a) to (c)

Quantification of Recurrence Plot

The recurrence quantification analysis (RQA), which measures the recurrence point density and the diagonal and vertical line structures of the RP, make it possible to identify and quantify the transitions between periodic, laminar, and chaotic states before PAF attacks; therefore, facilitating the possible prediction of PAF.[33] In this study, seven features were extracted from the RP of each five-min segment of each record to characterize different patterns. The related features are as follows:

Lmax: It measures the length of the longest diagonal line segment in the RP, excluding the main diagonal line. This is a very important recurrence variable because it is inversely related to the largest positive Lyapunov exponent.[12] Positive Lyapunov exponents gauge the rate at which trajectories diverge and are the hallmarks for dynamic chaos. Thus, the shorter the Lmax, the more chaotic (less stable) the signal will be.

Lmean: It measures the average length of the diagonal lines whose lengths exceed the certain threshold lmin:

where p(l) is the number of diagonal structures whose length is l. Lmean represents the average time that two segments of the trajectory are close to each other and can be interpreted as the mean prediction time.[30]

Entropy (ENTR): It is the Shannon entropy of the distribution of the length of diagonal segments:

where P(l) is the probability density of the diagonal structure whose length is l and it is defined as p(l)/sum(p(l)). ENTR reflects the complexity of the RP in respect of the diagonal lines.

Vmax: It measured the length of the longest line segment which is vertical.

Trapping Time (TT): It measures the average length of vertical structures:

where p(v) is the number of diagonal structures whose length is v. Here, vmin=2 is used. TT estimates the mean time that the system will abide at a specific state or how long the state will be trapped.

Recurrence trend (RT): It measures the variation rate of the REC away from the main diagonal which can reflect the drift and non-stationarity in a time series.

Firstly, the bottom right trigonal area of the RP is divided into K parts by K-1 equally spaced lines which parallel to the 45° diagonal, then the recurrence rate of the kth part REC is obtained and represents the mean of the sequence {REC,k=1,2,…,K}. RT is the linear regression slope of the sequence {REC, k=1,2,…,K} about {k=1,2,…,K} and it is calculated as follows:

Feature Dimension Reduction by LDA

In this step, we have used a feature dimension reduction technique to reduce the numbers of original features and to enhance the performance of the proposed algorithm.

LDA searches for those vectors in the underlying space that best discriminate among classes rather than those that best describe the data.[34] The goal of LDA is to seek a transformation matrix W that maximizes the ratio of the between-class scatter to the within-class scatter. Initially, we consider a within-class scatter matrix for the within-class scatter. A within-class scatter matrix S is defined as

where c is the number of classes, c is a set of data belonging to the ith class and m is the mean of the ith class. The within-class scatter matrix represents the degree of scatter within classes as a summation of covariance matrices of all classes. Next, we consider a between-class scatter matrix for between-class scatter. A between-class scatter matrix S is defined as

where m is the mean of all classes. The between-class scatter matrix represents the degree of scatter between classes as a covariance matrix of means of all classes. We seek a transformation matrix W that in some sense maximizes the ratio of the between-class scatter and the within-class scatter. The criterion function J(W) can be defined as

We can obtain the transformation matrix W as one that maximizes the criterion function J(W). Furthermore, given a number of independent features relative to which data is described, LDA creates a linear combination of these which yields the largest mean differences of the desired classes.[35] As a result, if there are c classes, the dimension of feature can be reduced to c-1 extremely.

In this paper, the numbers of original features have been reduced to three features by means of the LDA technique.

Classification

The last step is the classification of the HRV episodes by considering the extracted features. Different classification methods have been used for AF classification and prediction in the past.[935-40] In this work, the SVM classifier is used for classifying the episodes before PAF and far from PAF.

Given a training set (x,y),i = 1,2,…,l, where x ∈ R and, y ∈ {–1,1}, the traditional SVM algorithm is summarized as the following optimization problem:

subject to: y (w φ(x) + b) ≥ 1 – ξ , ξ > 0 ∀i.

where φ(x) is a nonlinear function that maps x into a higher dimensional space.[41] w, b, and ξ are the weight vector, bias, and slack variable, respectively. C is a constant and determined a priori. Searching for the optimal hyperplane in (10) is a quadratic programming problem, which can be solved by constructing a Lagrangian and transforming it into a dual maximization problem of the function Q(α), defined as follows:

where K(x, x) = φ(x) φ(x) is the Kernel function and, α = (α1, α2,…,α)is the vector of nonnegative Lagrange multipliers.

Assuming that the optimum values of the Lagrange multipliers are denoted as α (i = 1,2,…,l), it is then possible to determine the corresponding optimum value of the linear weight vector W and the optimal hyperplane as in (12) and (13), respectively:

The decision function can be written as:

In this work, the radial basis function is used as the kernel function and the parameters - kernel width σ and regularization constant C- were experimentally defined to achieve the best classification result.

RESULTS AND DISCUSSION

As we can see from Figures 1and 2, the RPs of the last 15- min of record far from PAF include many isolated recurrence points and the diagonal and vertical lines have short lengths. But in the record before PAF, as we are near to the end of record, we see some square-like patterns in the RPs and the length of diagonal and vertical lines become longer than previous minutes. In order to characterize these distinct patterns, we extracted the earlier mentioned features from the RPs of the fourth, fifth, and last five-min of each episode and considered the average value of each parameter to form the feature vector.

The boxplots of the features are presented in Figure 3 as we can see from this figure, the length of longest diagonal segment, Lmax, which inversely scales the largest positive Lyapunov exponent represents higher value in episodes before PAF. It can be interpreted that the episodes before PAF are more stable (less chaotic) than non-PAF episodes. The average length of the diagonal lines, Lmean, which reflects the average duration of a stable interaction in the signal, also exhibits higher value in the episodes before PAF. Entropy of the length of diagonal segments exhibits statistically significant increase in PAF episodes compared to non-PAF episodes. The increase of the value of entropy in PAF episodes is related to the more complexity of the deterministic structure of the RP before PAF events. Trapping time, which measures the average length of vertical structures, has higher value in PAF episodes compared to non-PAF ones. It means that the mean time that the system abides at a specific state in PAF episodes is statistically longer than non-PAF episodes. The length of the longest vertical segment, Vmax, also presents higher value in PAF episodes. Recurrence trend, which reflects the drift and non-stationarity in the HRV signal, also exhibits higher value in the episodes before PAF.

Figure 3

Boxplots of the features of episodes before PAF (PAF episodes) and episodes far from PAF (non-PAF episodes). (a) Lmean, (b) Lmax,(c) entropy, (d) trapping time, (e) Vmax, and (f) recurrence trend

To evaluate our proposed method, the features were extracted from AFPDB database. Afterward, the six original features are reduced to three new features by means of LDA. We used 50 ECG segments to train the SVM classifier and 56 segments to test the algorithm. The training and test sets belong to different subjects and there is no overlapping between the two sets.

To optimize the learning cost and the prediction performance, the SVM classifier parameters, kernel width σ and regularization constant C, must be chosen effectively. For this purpose, we have divided the training data into train and validation sets. The optimum values of the parameters are chosen such that the smallest error on validation dataset is achieved. The optimum values of parameters are achieved as 0.5 and 10 for σ and C, respectively.

To evaluate the performance of the proposed method, three measures are used as follows:

where TP, TN, FP, and FN stand for true positive, true negative, false positive, and false negative, respectively. If for example, a PAF episode is classified as PAF episode, then it is said that the episode is classified TP. On the other hand, if a non-PAF episode is classified as non-PAF episode, then it is said that the episode is classified TN. Any non-PAF episode which is classified as PAF episode by mistake will produce a FP; while any PAF episode which is classified as non-PAF episode by mistake will produce a FN result. The results of classification of the test data for both six original features and the reduced features are presented in Table 1. The obtained sensitivity, specificity, positive predictivity of the proposed method was 100%, 93.55%, and 100%, respectively.

Table 1

Comparison of the performances of original features and reduced features for PAF prediction

As seen from the results presented in Table 1, using LDA technique, not only reduces the number of the input features, but also increases the classification accuracy by selecting most discriminating features.

Several methods have been developed for automatic prediction of PAF using ECG signal in the past. Different automatic PAF prediction methods together with their reported results in terms of the commonly used measures of sensitivity and specificity are summarized in Table 2. All of these methods used AFPDB database to evaluate their approaches. The authors of[6363842] participated in Computers in Cardiology Challenge 2001 (PAF prediction challenge) and Zong et al.[6] could achieve the highest score. They used the number and timing of PACs as their main discriminator to predict PAF events and achieved a sensitivity of 79% for predicting the onset of PAF. Thong et al.[7] in 2004 developed an algorithm based on the analysis of isolated PACs not followed by a regular RR-interval in 30-min ECG segments. They calculated the difference in the number of PACs between two records of a subject and the record with the larger number of PACs considered as the one preceding PAF. They achieved a sensitivity of 89% and a specificity of 91% for predicting the onset of PAF. Chesnokov[9] in 2008 used complexity and spectral analysis of HRV, but his method couldn’t achieve good sensitivity for prediction of PAF events. Our proposed algorithm presents acceptable performances in terms of both sensitivity and specificity in identifying episodes preceding PAF. As seen in Table 2, the obtained results are better than the other existing approaches. Moreover, as in our method, there is no overlapping between the training and test datasets; it seems to be an efficient PAF predictor in real clinical applications.

Table 2

Summary of different automatic methods for PAF prediction and their reported results

CONCLUSION

In this paper, we have presented an effective HRV-based algorithm to predict the onset of PAF attacks. We have demonstrated that the structures of the RP contain useful information which can characterize distinct levels of chaoticity among episodes before PAF and episodes distant from PAF events. We have also shown that features based on the structure of RP including Lmax, Lmean, entropy, trapping time, Vmax, and recurrence trend can be used as good markers for predicting the onset of PAF attack. These extracted features are fed to a SVM classifier to identify episodes preceding PAF events. To evaluate the proposed method for predicting PAF, we have used the AFPDB database. The overall sensitivity, specificity, and positive predictivity of the proposed method in predicting the PAF events were 96.55%, 100%, and 100%, respectively. Comparing the performance of the proposed algorithm to those of the previously reported methods in the literature shows that the proposed algorithm is more effective than any of those methods.

BIOGRAPHIES

Maryam Mohebbi was born in Tehran, Iran, in 1981. She received the M.Sc. degree from the K. N. Toosi University of Technology, Tehran, Iran, and the Ph.D. degree from Tarbiat Modares University, Tehran, Iran, in biomedical engineering in 2007 and 2011, respectively. Dr. Mohebbi joined the Biomedical Engineering group at K. N. Toosi University of Technology in 2012.

Her research interests include biomedical signal processing, nonlinear analysis of ECG and HRV signals, modeling of cardiac activity, and pattern recognition.

Hassan Ghassemian was born in Iran in 1956. He received the B.S.E.E. degree from Tehran College of Telecommunication in 1980 and the M.S.E.E. and Ph.D. degree from Purdue University, West Lafayette, USA in 1984 and 1988 respectively. Since 1988, he has been with Tarbiat Modares University in Tehran, Iran, where he is a Professor of Computer and Electrical Engineering.

His research interests focus on Multi-Source Signal/Image Processing and Information Analysis and Remote Sensing. He is an IEEE Senior member.

Babak Mohammadzadeh Asl was born in Ahar, Iran, in 1982. He received the M.Sc. degree from the University of Tehran, Tehran, Iran, and the Ph.D. degree from Tarbiat Modares University, Tehran, Iran, in biomedical engineering in 2006 and 2010, respectively. Dr. Mohammadzadeh Asl joined the Biomedical Engineering group at Tarbiat Modares University in 2011.

His research interests include biomedical signal processing, array signal processing, and adaptive beamforming for medical ultrasound imaging.