Central tendency measure and wavelet transform combined in the non-invasive analysis of atrial fibrillation recordings

Paroxysmal AF database

Fifty Holter recordings of 30 seconds in length and two leads (II and V1) available in Physionet [19] were analyzed. The database included 26 non-terminating PAF episodes (group N), which were observed to continue in AF for, at least, one hour following the end of the excerpt, and 24 PAF episodes terminating immediately after the end of the extracted segment (group T). These signals were digitized at a sampling rate of 128 Hz and 16-bit resolution. Nonetheless, they were upsampled to 1024 Hz in order to allow better alignment for QRST complex subtraction, such as Bollmann et al. suggested [10]. This processing step is necessary to extract the AA from the surface ECG, see section Data preprocessing.

Persistent AF database

Sixty-three patients (20 men and 43 women, mean age 73.4 ± 9.0 years) with persistent AF lasting more than 30 days, undergoing ECV were followed during four weeks. A standard 12-lead ECG was acquired for each patient during the whole procedure and a segment of 30 seconds in length preceding the cardioversion was extracted from each recording for the analysis. All signals were digitized at a sampling rate of 1024 Hz and 16-bit resolution.

After the ECV, 22 patients (34.93%) maintained NSR during the first month. On the contrary, in 31 patients (49.20%), NSR duration was below 1 month and the remaining 10 (15.78%) relapsed to AF immediately after ECV. These 41 patients constituted the group of AF recurrence. All patients were in drug treatment with amiodarone. The median arrhythmia duration was 10.58 months (range 1–47.22) and echocardiography demonstrated a mean left atrium diameter (LAD) of 45.82 ± 6.93 mm. 20.63% of the patients presented underlying heart disease. No statistically significant differences were found in the aforementioned clinical parameters between the patients who maintained NSR and relapsed to AF.

Data preprocessing

In both databases, lead _{V 1}was chosen for the analysis because previous works have shown that AA is dominant in this lead [20]. This signal was preprocessed using forward/backward highpass filtering with 0.5 Hz cut-off frequency to remove baseline wander, next, lowpass filtering with 70 Hz cut-off frequency was applied to reduce high frequency noise and, finally, notch filtering at 50 Hz was applied to remove powerline interference [21]. On the other hand, reliable analysis of the AA from surface ECG recordings requires that ventricular activity has first been cancelled [10]. Although a variety of different techniques exist for this purpose, an adaptive singular value QRST cancellation template was applied [22].

Wavelet Transform

where a b∈ℜ a≠0 are the scale and translation parameters, respectively, and t is the time. As a increases, the wavelet becomes narrower. Thus, one have a unique analytic pattern and its replications at different scales and with variable time localization.

The decomposition results in wavelet coefficients C, which depend on scale and position. In fact, a vector of wavelet coefficients _{C m} is obtained for each analyzed discrete scale m. The information stored in the wavelet coefficients vectors is not repeated elsewhere and allows the complete regeneration of the original signal without redundancy, because the used discretization of the mother wavelet employs orthonormal basis functions [23].

Central Tendency Measure

Chaotic equations are sometimes used to generate graphs. Thus, the graph x(n + 2)−x(n + 1) versus x(n + 1)−x(n) produces a scatter plot of first differences of the data, x(n) being the value of a time series at time n. This plot, centered around the origin, gives a graphical representation of the variability degree in the time series and is useful in modeling biological systems, such as hemodynamics and heart rate variability [14]. With this approach, rather than defining a time series as chaotic or not chaotic, the degree of variability or chaos is evaluated.

In the present work, f waves organization was estimated by computing the CTM from the first differences scatter plot of the wavelet coefficients vector for the scale containing the typical AF frequency range. For this purpose, a seven-level wavelet decomposition was applied to the AA signal. The seventh discrete scale resulted in a good match to the AA frequency band of interest [24], i.e. 4–8 Hz. In fact, the dominant atrial frequency (DAF), i.e., the highest amplitude frequency, was within this range for all the analyzed patients. Furthermore, previous works have proved that values for this frequency below 4 Hz or above 8Hz are unusual [25].

A good matching between wavelet scale and AF frequency range is relevant to preserve the essential behavior of the arrhythmia. In this respect, the DAF has provided to be a concomitant indicator of atrial refractroriness [26], thus revealing very clinically interesting information about AF progression under different therapies [16]. Remark that for sampling rates different from 1kHz, the seventh wavelet scale could not match properly to the typical AF frequency range and, consequently, the proposed method performance could be notably reduced. Thus, resampling up to 1kHz of new input ECG recordings has to be considered as a previous step before their processing under the proposed method.

Optimal parameters selection

The proposed method performance depends on the chosen wavelet function and the radius ρ selected for computation of CTM. However, there are no established rules for the choice of their optimal values [2]. To this respect, a cautious and still exploratory approach is to test different wavelet families and then to compare their efficiency in the specific problem [27]. Unfortunately, on each electrocardiographic application where the WT has been used, a different wavelet family was chosen [4]. In this study, several orthogonal wavelet families were tested, because only in an orthogonal basis any signal can be uniquely decomposed and the decomposition can be inverted without loosing information [23]. Thus, all the different functions from Haar, Daubechies, Coiflet, Biorthogonal, Reverse Biorthogonal and Symlet wavelet families were tested.

On the other hand, an approach similar to the developed by Hornero et al[28] was used to select the optimal ρ value. Thus, the CTM was computed for radius values of 0.1, 0.2, …, 10 times the standard deviation of the analyzed data, i.e., the wavelet coefficients vector associated with the seventh discrete scale. Normalizing ρ in this way provides translation and scale invariance, in the sense that CTM remains unchanged under uniform process magnification, reduction or constant shift to higher o lower values. To this respect, note that the wavelet vector interesting information is carried by the coefficients temporal variations, rather than its absolute amplitude variations. In fact, the amplitude is only dependent on the similarity between the AA under study and the selected mother wavelet. Hence, amplitude gives no useful information about f waves variability.

Statistical analysis

As the considered number of episodes for both analyzed databases was not notably large, the proposed method was validated by a resubstitution approach, i.e., it was learned from all the available data and then tested on the same set of data. The CTM value, for a specific ρ value and wavelet function, providing maximum discrimination between groups, that is, the optimum threshold, was obtained by means of a receiver operating characteristic (ROC) curve. The ROC curve is a graphical representation of the trade-offs between sensitivity and specificity when discrimination threshold is varied. For PAF patients, sensitivity (i.e. the true positive rate) was considered as the proportion of non-terminating PAF episodes correctly discerned, whereas specificity (i.e. the true negative rate) represented the percentage of terminating episodes properly identified. Similarly, for the ECV outcome analysis, sensitivity was the proportion of patients relapsing to AF appropriately classified and specificity was the percentage of patients resulting in NSR after ECV accurately predicted. The total number of PAF patients and ECVs precisely classified was considered as the diagnostic accuracy corresponding to each prediction. The CTM value providing the highest accuracy was selected as optimum threshold. The area under the ROC curve (AROC) was also computed. AROC is a single number which represents a summary of performance. For a perfect test the area is 1, while an AROC of 0.5 represents a worthless test.

To verify representativeness of the used databases and generality of the results obtained by the resubstitution approach, a leave-one-out cross-validation (LOOCV) scheme was also used. In this procedure, the method was trained and tested a number of times equal to the number of patients in the database. In each iteration, all the data, except for a single observation, were used for training and the model was tested on that single observation. In addition, significant differences between terminating and non-terminating PAF episodes and between patients who resulted in NSR and relapsed to AF were evaluated making use of a Student’s t-test. All the groups had a normal and homoscedastic distribution as the Shapiro–Wilk and Levene tests proved, respectively. A two-tailed value of statistical significance p<0.05 was considered statistically significant.

Optimal parameters selection

After optimal wavelet function selection, the proposed method dependence on the radius ρ was investigated following the process described in section Optimal parameters selection. For each tested ρ value, statistical differences between terminating and non-terminating PAF episodes and between patients who relapsed to AF and maintained NSR after ECV were evaluated by using Student’s t-test. In both cases, statistically significant differences were found for ρ values between 2.8 and 5.6. Nonetheless, the lowest p value was noticed for ρ=3.3 and ρ=4 in PAF termination and ECV outcome predictions, respectively. Thereby, the CTM was computed with these ρ values in each analyzed AF scenario.

Classification performance

For PAF termination prediction, CTM values computed with previously described optimal parameters provided sensitivity, specificity and accuracy of 100% (26 out of 26), 91.67% (22 out of 24) and 96% (47 out of 50). In addition, 94% (46 out of 50) of cross-validated grouped cases were correctly identified. The ROC curve provided 0.888 as optimum discrimination threshold between terminating and non-terminating PAF episodes with an AROC of 0.974, see Figure 1(a). Additionally, Figure 1(b) shows that the non-terminating PAF episodes presented higher CTM values (0.963 ± 0.036) than the terminating ones (0.855 ± 0.040). As an example, Figure 2 shows a 10 second-length ECG interval, its extracted AA signal, the wavelet coefficients vector corresponding with the seventh discrete scale and its first differences scatter plot for a typical terminating PAF episode and other non-terminating one. Faster temporal variations within the wavelet coefficients vector and, hence, higher dispersion in scatter plot and lower CTM value can be appreciated for the terminating AF episode.

Regarding ECV outcome prediction, the ROC curve provided 0.823 as optimum CTM discrimination threshold, in which 82.93% (34 out of 41) sensitivity and 90.91% (20 out of 22) specificity were obtained, see Figure 3(a). Hence, the ECV outcome in 54 out of 63 patients (85.71%) was correctly predicted, with 82.54% (52 out of 63) of cross-validated grouped cases appropriately identified. An AROC of 0.940 was achieved. Moreover, the patients relapsing to AF presented lower CTM values (0.761 ± 0.054) than those resulting in NSR after one month (0.880 ± 0.050), such as Figure 3(b) shows. As for PAF termination prediction, Figure 4 presents typical ECG intervals, extracted AA signals, wavelet coefficients and scatter plots for a patient who maintained NSR and other relapsing to AF during the first month post-cardioversion of persistent AF. In this case, a higher variation rate within the wavelet coefficients, associated with a lower CTM value, was observed for patients who relapsed to AF.

Other interesting result from Figure 3(b) is that all the patients relapsing to AF immediately after ECV were appropriately identified with a CTM threshold of 0.749. Indeed, these patients presented CTM values (0.727 ± 0.021) lower than the remaining ones who relapsing to AF during the first month post-cardioversion (0.771 ± 0.057), the statistical significance being lower than 0.001. Nonetheless, note that the false positive rate was 32.26% (10 out of 31).

Comparison with previous works

Non-linear analysis metrics are a valuable tool in the assessment of physiological time series because hidden information related to underlying mechanisms can be obtained [31, 32]. To this respect, it has become evident through a number of studies that the cardiac system is non-linear on its behavior [33]. Furthermore, regarding the employment of non-linear metrics, Sample Entropy (SampEn) has been successfully applied to the surface ECG in AF. In this case, SampEn has been computed over the main atrial wave (MAW). The MAW is the fundamental waveform associated to the AA signal and has served as the starting point to estimate non-invasively AF organization. Through the MAW-SampEn strategy, it has been possible to achieve a diagnostic ability of 90% in the prediction of PAF termination [34] and approximately 80% in the prediction of ECV outcome [35]. Similarly, f waves regularity has also been assessed through the application of SampEn to the wavelet domain of the AA signal. In this case, the wavelet-SampEn approach reached a discriminant capability of 90 and 82% for PAF termination and ECV outcome predictions, respectively [29, 30]. Nonetheless, SampEn might be notably sensitive to the presence of noise [34, 36] and spikes [37]. Thus, given that QRST residua are often present in the extracted AA [20, 22] and the predictive ability of Wavelet-SampEn could be improved, the analysis of alternative non-linear methods is desirable to develop top performance single classificators of AF organization-related events.

Regarding the frequency domain, a variety of works analyzed the DAF of the AA signal, providing diagnostic abilities in PAF termination prediction between 86% and 90%, depending on the used method for its computation [40, 41]. Holmqvist et. al. [44] evaluated a parameter obtained from time-frequency analysis of the AA signals, such as harmonic decay [24], but a low number of patients relapsing to AF after ECV were correctly identified (47%). Finally, in [46], ventricular rhythm was analyzed using three-dimensional RR intervals plots, quantifying clustering of RR intervals, however, only 50% of patients who relapsed to AF during the first 4 weeks following ECV were correctly discerned.

Nevertheless, single parameters combinations and advanced classification tools can be found in the literature to predict AF behavior. Thus, regarding PAF termination prediction, DAF has been combined with other parameters to improve its diagnostic accuracy. To this respect, Petrutiu et al[42] studied the AA peak frequency power evolution within the last two seconds before spontaneous PAF termination, reaching thus an accuracy of 93.33%. The same classification result was reported by combining SampEn and WT [29]. In that work, SampEn was used to assess f waves regularity from the time and wavelet domains of the AA signal, providing two different and independent classifications, which were combined as a function of the DAF. A slightly higher result (accuracy of 96.67%) was reached by Sun & Wang [39] making use of a multilayer perceptron neural network to combine 11 features extracted from the ECG recurrence plot quantification. With regard to ECV outcome prediction, Watson et al[43] examined a variety of wavelet transform-based statistical markers, which obtained a sensitivity of 88% and specificity of 100% by means of a non-parametric classification system. Zohar et al[45] developed a non-deterministic model with several parameters as inputs for predicting NSR maintenance after ECV, providing a diagnostic accuracy of 84%. As for PAF termination prediction, the AA signal organization estimation both in time and wavelet domains through SampEn provided an sensitivity of 95% and specificity of 93% [30]. However, it is worth noting that in this study only patients undergoing the first attempt of ECV and, therefore, with AA signals notably organized [47], were analyzed. Finally, recent works have reported that the combination of f waves amplitude with the DAF [38] and SampEn [35] computed from the MAW reached a discriminant ability of 86 and 90%, respectively.

Bearing this context in mind, it can be noticed that only complex combinations of single predictors can improve the CTM classification result. Thereby, the proposed metric can be considered as a promising single estimator of PAF termination and ECV outcome prediction, with the additional advantage of a simpler implementation to work in real-time. To this respect, DWT can be computed efficiently with a pyramid filter bank algorithm [23], allowing thus its implementation in real-time [48, 49]. In contrast, the classical algorithm proposed in the literature for SampEn computation requires a high execution time, which is not fast enough for online applications [50]. Although faster alternatives for SampEn computation have been recently proposed [50, 51], their ability to be implemented in real-time environments has not been proved yet. In addition, accuracy of these new algorithms has not been validated by comparison with the classical SampEn definition. On the other hand, it has to be remarked that the complex combination of multiple parameters or the use of advanced classification techniques, such as in [39] or in [43], makes difficult the clinical interpretation of the results. In this sense, possible clinical meaning of each parameter is blurred within the classification approach.

Interpretation of results

Taking a specific wavelet decomposition scale, the coefficients vector contains the similarity evolution through time between the analyzed signal and the scaled mother wavelet. A low variability in this time series implies an invariable waveform, regardless of its concrete shape, along the studied time period. In contrast, high variability implies variable waveforms that may, eventually, evolve to a more organized pattern. This fact justifies the results obtained in the prediction of spontaneous PAF termination, where terminating episodes presented wavelet vectors with higher variability than non-terminating ones, see Figure 1. The AA evolution from disorganized f waves to organized P waves that takes place in AF recordings prior to its termination may cause this behavior in the wavelet coefficients vector [20].

Regarding ECV result analysis, the presence of more structured f waves in organized AA signals [15, 20] could justify the obtained results, which show that patients who relapsed to AF presented wavelet coefficients vectors with higher variability than those who remained in NSR, see Figure 3. In fact, patients relapsing to AF immediately after ECV presented, in average, lower CTM values than the remaining ones. These findings, suggesting more organized AA signals in effective cardioversions one month after the procedure, agree with observations obtained from previous works, such as: (i) the higher the AA organization, the higher the success rate in AF cardioversion [17, 52], (ii) the higher the AA organization, the lower the energy required for successful cardioversion [18] and (iii) PAF requires less energy for cardioversion than persistent AF [53]. These observations highlight the fact that, when a higher number of reentries are wandering throughout the atrial tissue, a lower probability of successful ECV is obtained. One possible explanation could be that a low degree of AA organization might result in an increased mass of atrial myocardium that is not fully excitable [52].

Study limitations

The analysis was developed with a limited group of patients. Hence, a larger sample allowing a more rigorous statistical study would be required to provide improved confidence in the robustness of the developed approach. To this respect, wider databases containing non-terminating and spontaneously terminating PAF episodes after different time epochs (ten minutes, half an hour, an hour, ten hours, etc.) and patients who resulted in NSR and relapsed to AF after 3, 6, and 12 months following ECV would be necessary. In addition, the availability of longer PAF recordings would allow to address the interesting question about the time in advance with which spontaneous PAF termination could be predicted. On the other hand, the persistent AF database only included suitable ECV patients following the standard clinical criteria; therefore, it is unknown how the proposed algorithm will behave in patients with adverse clinical predictors, like atrial dilatation, etc., which, by default, are excluded from ECV procedures. Finally, only lead _{V 1} was analyzed rejecting the possible information contained in the remaining leads. However, for this type of studies, lead _{V 1}seems to be the most suitable lead because significant correlations between the atrial frequency [54] and SampEn [15] obtained from this lead and those obtained from atrial electrograms have been observed.