Machine Learning Approach to Predict Ventricular Fibrillation Based on QRS Complex Shape

Early prediction of the occurrence of ventricular tachyarrhythmia (VTA) has a potential to save patients’ lives. VTA includes ventricular tachycardia (VT) and ventricular fibrillation (VF). Several studies have achieved promising performances in predicting VT and VF using traditional heart rate variability (HRV) features. However, as VTA is a life-threatening heart condition, its prediction performance requires further improvement. To improve the performance of predicting VF, we used the QRS complex shape features, and traditional HRV features were also used for comparison. We extracted features from 120-s-long HRV and electrocardiogram (ECG) signals (QRS complex signed area and R-peak amplitude) to predict the VF onset 30 s before its occurrence. Two artificial neural network (ANN) classifiers were trained and tested with two feature sets derived from HRV and the QRS complex shape based on a 10-fold cross-validation. The prediction accuracy estimated using 11 HRV features was 72%, while that estimated using four QRS complex shape features yielded a high prediction accuracy of 98.6%. The QRS complex shape could play a significant role in performance improvement of predicting the occurrence of VF. Thus, the results of our study can be considered by the researchers who are developing an application for an implantable cardiac defibrillator (ICD) when to begin ventricular defibrillation.


INTRODUCTION
Ventricular tachyarrhythmia (VTA) causes a rapid heart rate and eventual death in the absence of immediate medical intervention (Lee et al., 2016). As the majority of sudden cardiac deaths (SCD) occur because of VTA (Lee et al., 2016), early prediction of VTA is important to save patients' lives. VTA contains different types of arrhythmias, such as ventricular tachycardia (VT) and ventricular fibrillation (VF). Because measuring and analyzing electrocardiogram (ECG) signals is an efficient way to identify electrical conduction malfunctions in the heart, such as arrhythmias, previous studies have attempted to predict the occurrence of VT, VF, or both by investigating electrocardiography (ECG) (Riasi and Mohebbi, 2013;Cappiello et al., 2015;Melillo et al., 2015).
Various methods have been introduced to predict VTA (VF, VT, or both), such as by assessing QRS (Q, R, and S waves in ECG) duration, T wave alternans, left ventricular impairment, QT (from the start of Q wave to the end of T wave in ECG) dispersion, and heart rate variability (HRV) (Lane et al., 2005). Among these, HRV is the most commonly employed signal that provides features for isolating arrhythmia from the normal HRV (Reed et al., 2005). HRV is a measure that indicates time variation in consecutive heartbeats, it is also denoted as RR (Billman et al., 2015). HRV has been analyzed to quantify its features using three analysis methods: time domain, frequency domain, and Poincare non-linear analyses (Bilgin et al., 2009;Joo et al., 2012;Lee et al., 2016).
Previous studies mainly used the above-mentioned three analysis methods to predict VT, VF, or both using HRV. Bilgin et al. (2009) performed a sub-band frequency analysis on HRV data and presented its feasibility on VTA prediction as compared to the traditional frequency analysis using two base-bands; low frequency (0.04-0.15 Hz) and high frequency . To obtain the sub-bands, they used a wavelet packet transform (WPT) and evaluated the sub-bands using a multilayer perceptron (MLP) neural network. Elias et al. used various features extracted in time and frequency domain (Ebrahimzadeh and Pooyan, 2011) to detect SCD early in patients with sustained VTA. They used an MLP neural network and k-nearest neighbors (KNN) to classify the healthy subjects and those prone to SCD, and principal component analysis (PCA) to reduce the feature dimensions. Recently, Elias et al. used both time-frequency and Poincare non-linear analyses to extract HRV features (Ebrahimzadeh et al., 2014). To evaluate the performance of their methods in the prediction of SCD in patients with sustained VTA, the features were extracted from different segments of HRV signals at successive 1-min intervals (i.e., the first, second, third, and fourth minute before the event). MLP and KNN were again used to classify healthy and VTA (Ebrahimzadeh et al., 2014). Joo et al. (2012) predicted VT and VF 10 s before their occurrences using HRV features. They applied an artificial neural network (ANN) to deal with the complexities of the features extracted from the HRV (Joo et al., 2012). Recently, Lee et al. (2016) used the respiratory rate variability (RRV) features to improve the accuracy of VTA predictions using ANN. The performance of their predictions using RRV features outperformed HRV features.
Several studies showed promising performances of predicting VTA using HRV, RRV, other features, and their combination (Lee et al., 2016). However, they did not consider QRS complex features. The QRS complex represents the electrical activation of ventricles (Zhang et al., 1997) from which important features can be extracted. If more emphasis is placed on feature extraction from QRS complexes, performance of predicting VTA could be improved significantly.
The shape of the QRS complex provides abundant information about the pattern of the electrical propagation through ventricular tissue (Zhang et al., 1997). In clinical applications, the feature extraction and analysis of QRS complexes can predict ventricular arrhythmia, e.g., VF (Bassareo and Mercuro, 2013). Therefore, we hypothesize that the features from QRS complexes could be used to predict VTA in advance. The aim of this study was to investigate the feasibility of using the features extracted from QRS complexes for the early prediction of VTA (i.e., VF), as compared to traditional HRV features. To this end, we extracted two features such as QRS singed area and R-peak amplitude and investigated the prediction performance of VF using ANN, anticipating increased prediction performance using the features from the QRS complex. Also, four alternative machine learning algorithms showed a similar trend as ANN with high prediction performance using QRS shape features.

Dataset
We used the following freely available databases in PhysioNet (RRID:SCR_007345) (Goldberger et al., 2000): Creighton University (CU) ventricular tachyarrhythmia (CUDB) (Nolle et al., 1986), normal datasets from paroxysmal atrial fibrillation (PAF) prediction challenge database (PAFDB) (Moody et al., 2001) and the MIT-BIH normal sinus rhythm database (NSRDB) (Goldberger et al., 2000). The sampling frequencies were 250 Hz for the CUDB and 128 Hz the other two databases. Although there were 35 recordings in the CUDB, seven recordings were not considered because some contained only VT events (not VF events) and others were shorter than the required data length (>150 s). Thus, a total of 27 recordings were used for data analysis. The control group consisted of 28 subjects (22 subjects who did not have fibrillation events from the PAFDB and 6 subjects from the NSRDB).

Preprocessing
The databases in PhysioNet (PhysioNet, RRID:SCR_007345) (Goldberger et al., 2000) contained raw ECG signals and their corresponding HRV signals. The R-peak to R-peak intervals (RR) were produced by reading the annotation files in PhysioNet (PhysioNet, RRID:SCR_007345) (Goldberger et al., 2000) that were annotated by Cardiologists. Figures 1A,B show examples of ECG and HRV signals before VF occurred, respectively. We divided the signal into two parts: required and forecast time. The required time represents the time period used for feature extraction between 150 and 30 s before the VF onset time. The forecast time is the time period between 30 and 0 s before VF onset. Using the required time data, we could predict the occurrence of VF before the forecast time.

Feature Extraction
Features were extracted from 27 VF and 28 control datasets. The descriptions of the features used in this study are listed in Table 1 (Lee et al., 2016). These features consist of 11 FIGURE 1 | (A) ECG and (B) HRV signals of a VF subject that include two pre-VF regions. The signal 30 s before the VF onset is called the forecast time and the signal 120 s before the forecast time is called as required time. RR represents R-peak to R-peak interval in milliseconds (ms).
HRV features (4 features in the time domain, 4 features in the frequency domain, 3 features using Poincare non-linear analysis) and 4 QRS complex features (in the time domain). The desired features for the investigation were extracted from the required time between 150 and 30 s before VF onset. The QRS complex signed area and R-peak amplitude were computed using a function known as ECG-derived respiratory (EDR) from the PhysioNet (PhysioNet, RRID:SCR_007345) MATLAB toolbox (Moody et al., 1985). The method in EDR uses the effect of modulation of the R-peak amplitude, which is evaluated by processing the signed area under QRS complex in the ECG signal (Mazzanti et al., 2003). We computed the QRS complex signed area using the weighted sum of the samples between two boundaries: the interval from the PQ junction (the point where P wave and Q wave meet in ECG) to J-point (beginning of ST segment) (Moody et al., 1985). The R-peak amplitude was also computed by selecting a sample that had maximum value between the boundaries. The RR (in the equations) represents the R-peak to R-peak interval. In this end, the following 4 QRS features were used for early prediction of VF: mean and standard deviation of the QRS complex signed area and the R-peak amplitude.

HRV Features
All HRV features were computed from successive RR intervals.

Time Domain Features
Four HRV features were computed in this category: (1) mean RR intervals [Mean NN (RR)], (2) standard deviation of NN (RR) intervals (SDNN), (3) square root of mean squared difference of successive NN (RR) intervals (RMSSD), and (4) the proportion of interval differences of successive NN (RR) intervals greater than 50 ms by the total number of NN (RR) intervals (pNN50), defined as follows:

Frequency domain features
We considered three frequency bands, such as the very low frequency (VLF) band (0-0.04 Hz), low frequency (LF) band (0.04-0.15 Hz), high frequency (HF) band (0.15-0.4 Hz), and the ratio of LF and HF. We computed the power spectrum density RPampSD µV Standard deviation of the R-peak amplitude FIGURE 2 | The architecture of our ANN using 11 HRV features. The input features to the ANN. Mean NN: mean normal R-peak to normal R-peak interval, SDNN: standard deviation of NN, RMSSD: square root of mean squared difference of successive NN, pNN50: proportion of interval differences of successive NN intervein greater than 50 ms, VLF: very low frequency, LF: low frequency, HF: high frequency, SD1: standard deviation of points perpendicular to the axis of line of identity, SD2: standard deviation of points along the axis of identity, and the ratio of SD1 and SD2.
(PSD) of the bands using Welch's periodogram with a Hanning window (window size: 256 points with an overlap of 50%).

Poincare non-linear features
The Poincare non-linear features were dispersion of points perpendicular and points along the axis of the line-of-identity. The standard deviation of the successive RR intervals scaled by 1/v2 (SD1) and the standard deviation of points along the axis of line-of-identity (SD2) were both calculated using (5) and (6). We considered the ratio of SD1 and SD2 as well.

QRS Complex Features
The QRS complex shape includes Q, R, and S waves from which the signed areas and the R-peak were calculated. The mean for QRS complex signed area and R-peak amplitude of the ECG were calculated using (7) and (8), and their standard deviations were calculated using (9) and (10).

Prediction Algorithms
The architecture of our ANN was a fully connected network structure consisting of three layers: an input layer with nodes representing input variables to the problem, a hidden layer containing nodes to help capture the nonlinearity of the input data, and an output layer with a node representing the dependent variable (Figures 2, 3) (Lippmann, 1989;Basheer and Hajmeer, 2000). The hidden layer consisted of six neurons which were selected by trial and error (Figures 2, 3) with rectified linear unit (RELU) (Glorot et al., 2011) activation functions, and the output layer used a sigmoid activation function (Narayan, 1997). Activation functions decide which neurons should be activated or deactivated (Leshno et al., 1993). We implemented two ANN models with two different input parameters: 11 HRV features (Figure 2) and 4 QRS shape features (Figure 3). We randomly shuffled the input features, and then used StandardScaler function from sklearn preprocessing library to Frontiers in Physiology | www.frontiersin.org standardize the input features. The features were standardized by removing the mean and scaling to unit variance. The standard score Z of a feature x is calculated as: z = (x−µ)/s, where µ is the mean and s is the standard deviation of input features of all datasets.
Hence, the governing equations for RELU (12) and sigmoid (15) activation functions are provided as following: f (X j ) = 0 for X j ≤ 0 X j for X j ≥ 0 (12) where x i is an input feature to the hidden layer, w ij is a weight of the connection between ith input and jth neuron in the hidden layer, X j is a weighted sum of the dot products of the input x j and weight w ij , and n is number of input features. The output from jth neuron of the hidden layer is y j (13), which is computed by applying RELU activation function to X j and adding bias b j . Also, the output y k of the output layer is computed by applying sigmoid activation function to X k , which is the weighted sum of the dot products of y j and a weight of the connection between jth neuron in the hidden layer and kth neuron in the output layer w jk , and adding bias b k . The m is number of neurons in the hidden layer. The ANN was trained using Adam optimizer algorithm which updates the weights (Kingma and Ba, 2014). Adam is an adaptive learning rate optimization algorithm that has been designed for training neural network. Because we implemented a binary classification model, we used binary cross-entropy as loss function to measure the divergence between two probabilities distribution.
Furthermore, we implemented four more machine learning algorithms such as: support vector machine (SVM), KNN, random forest (RF), and Gaussian Naive Bayes (NB) classifiers. These machine learning algorithms were implemented using sklearn library in python3.
All algorithms were evaluated 10 times with a 10-fold cross validation, to avoid overfitting. In the 10-fold cross validation (Weiss and Kulikowski, 1991), the dataset was randomly divided into approximately 10 groups. One group was treated as the testing dataset, and the remaining 9 groups were used for training. The cross-validation was repeated 10 times.
Finally, the prediction accuracies were estimated by calculating the means and standard deviations of each model. To observe the statistical differences between HRV and QRS shape accuracies, we performed two tailed t-test for each model. Also, to check the statistical differences among the accuracies for all the machine learning algorithms, we computed one-way repeated-measures ANOVA with Tukey post hoc analysis for multiple comparisons. The flowcharts of the methods used in this study are shown in Figure 4, where Figures 4A,B are based on 11 HRV and 4 QRS complex features, respectively. Table 2 shows the comparison of the means and standard deviations of the HRV and QRS complex shape features between the control and VF dataset (see Supplementary Materials for more illustration). Nine of the 15 features, SDNN, RMSSD, pNN50, LF, SD1, SD2, QRSaSD, QRsampM, and QRSampSD, show statistically significant differences (two tailed t-test, p < 0.05). Table 3 summarizes the performance of two ANNs with different feature types; HRV vs. QRS. Eleven HRV features achieved 72% prediction accuracy. The sensitivity and specificity were 65.68% and 98.44%, respectively. When using 4 features extracted from the QRS complex singed area and the R-peak amplitude, the prediction performance improved dramatically. The accuracy, sensitivity, and specificity were 98.6, 98.4, and 99.04%, respectively. The result shows that the QRS complex shape features extracted from the ECG could have an impact in predicting VF before its occurrence in terms of its prediction performance. Table 4 presents the average computational times required for training and testing ANN using HRV and QRS complex shape features. The computational times needed for training and testing ANN using 11 HRV parameters were 1545 and 0.72 ms, respectively. Similarly, the computational time needed for training and testing ANN using 4 QRS complex shape parameters were 1505 and 0.7 ms, respectively. Although the number of input features are different, there was no significant difference between the two ANN models in terms of computational time.  Note that the training time was estimated while an ANN model was constructed, and the testing time was estimated while one sample produced a prediction result. The training and testing times were computed for each cross-validation step, and they were averaged. Figure 5 presents the receiver operating characteristic (ROC) curve for the three models. The ANN with 11 HRV features has the lowest area under curve (AUC) value (0.71). When using 4 QRS complex shape features, the AUC reached 0.99. Figure 6 shows the means and standard deviations of the accuracies evaluated 10 times using a 10-fold cross validation for all machine algorithms we considered. The performances of the QRS shape features statistically outperform those of the HRV features for all machine algorithms (two tailed t-test,

DISCUSSION
In this study, we showed that the performance for predicting VF can be improved by using features extracted from the QRS complex shape (QRS complex signed area and R-peak amplitude). A maximum prediction accuracy of 98.6% was obtained using ANN when only 4 features of the QRS complex signed area and R-peak amplitude were used. However, using HRV features presented a significantly low performance with a prediction accuracy of 72%. The ROC curve in Figure 5 also showed the same trend. From the result in Table 3, we demonstrated that QRS shape features can predict VF before its occurrence more accurately than the traditional HRV features. As depicted in Figure 6, the prediction performance obtained using 4 QRS shape features was statistically higher than that was obtained using 11 HRV features for all algorithms. Also, ANN statistically outperformed three algorithms, such as KNN, RF, and NB, when using 4 QRS features, and its performance was comparable with SVM. However, ANN needed computational time more than the other algorithms to train the model, but not for testing time; all algorithms, including ANN, need less than 1 ms for one testing (Supplementary Table 2).
The QRS complex of an ECG contains information about the ventricular depolarization process. The time at which the QRS complex is generated is the time required to complete ventricular depolarization. The QRS amplitude is proportional to the energy consumed for ventricular depolarization. Therefore, if the electrical properties of the ventricular tissue remain unchanged, the QRS shape does not change within a short period of time (several seconds to several hours). How could we predict the VTA shortly before using QRS shape information? To our best knowledge, the reason is that when the ventricles begin to make reentrant waves due to ectopic focus or other reasons, the cardiac electrical wave pattern begins to change (which can affect the QRS shape). In addition, when electrical waves are different, the way in which the ventricles contract is also different, which changes the location of the ventricular tissue that is the source of the ECG. Changes in the location of ventricular tissue will affect the ECGs that reflect this. Therefore, QRS complex shape represents the depolarization (activation) of the ventricle muscles (Zhang et al., 1997) from which abnormalities in the electrical activation features can be extracted for early prediction of VF. The findings highlight the importance of the QRS complex shape features for predicting VF.
The means and standard deviations presented in Table 2 show the comparison between the features of VF and control datasets. Nine features have statistically significant differences between the VF and control groups (p < 0.05). We could use these features to compare the VF and control groups. However, the prediction of VF is not achievable by mere comparison of some parameters but by classification of complex patterns of every feature based on machine learning techniques.
Previous studies dealt with features extracted from HRV and predicted VF, VT, or both with promising performance. Elias et al. showed a performance with 99.73% accuracy for features extracted from HRV 1 min just before the occurrence of SCD early in patients with sustained VTA (Ebrahimzadeh et al., 2014). However, they did not consider the forecast time which is a time period before the occurrence of the VTA. Joo et al. (2012) predicted VT and VF 10 s before the FIGURE 6 | Means and standard deviations of the prediction accuracies of each algorithm. Single asterisk ( * ) indicates a statistically significant difference between prediction accuracies using HRV and QRS shape features (QRS > HRV, p < 0.001), and double asterisk ( * * ) between the prediction accuracies of different algorithms (ANN = SVM > KNN, RF, and NB, p < 0.001). HRV: heart rate variability with 11 parameters, QRS shape: QRS signed area and R-peak amplitude with 4 parameters, ANN: artificial neural network, SVM: support vector machine, KNN: k-nearest neighbors, RF: random forest, NB: Gaussian Naïve Bayes. Results show that features extracted from HRV contain important information for predicting the occurrence of VF several minutes in advance. However, Lee et al. (2016) revealed that the performance using only HRV features can be improved by adding RRV features. We found that only the QRS complex shape or that combined with HRV can improve the performance of predicting VF. In our study, we used 2-min-long signal to predict VF 30 s before its occurrence. The signals we used for the analysis (required time) and the prediction time gap (forecast time) were short. However, our study showed that the features extracted from QRS complex morphology (shape) could have effects for predicting VF.
We compared the performance obtained using a combination of HRV and QRS shape features with that obtained using only QRS shape features, but little improvement in prediction accuracy (only ∼1%) was found for the combination features. This indicates that using QRS shape features solely would be an efficient way to predict VF. Therefore, we decided to not include the result for the combination features in our study.
Our algorithm could be installed in patients' implantable cardiac defibrillator (ICD) for real-time VF prediction as an additional functionality to VF detection. Predicting the occurrence of VF hours in advance would be more useful, however, the datasets used for this paper limited to 120 s data window and predict VF 30 s before its occurrence. Au-Yeung et al. (2018) showed that a correct prediction could be made when the ventricular arrhythmia occurs nearer. Thus, our prediction accuracy of 98.6% was higher than that of Bayasi et al. (2016) who predicted the occurrence of VF 3 h prior to the onset with an accuracy of 86%, and Lee et al. (2016) who predicted the onset 1 h prior with an accuracy of 85.3%.
According to the Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, a recording of approximately 1 min is needed to calculate the HF component, and at least 2 min, to calculate the LF component (Task Force of the European Society of Cardiology the North American Society of Pacing Electrophysiology, 1996). However, the VLF calculated from short-term recordings (<5 min) is an uncertain measure (Task Force of the European Society of Cardiology the North American Society of Pacing Electrophysiology, 1996). The VF datasets in this study were very short in length, therefore, we had to use 120 s data window for extracting features 30 s before VF occurs.
The limitation of our study was the small dataset (55 recordings) and the short length of the signals before the VF occurred. To implement a study for clinical purposes, our ANN model must be trained using more datasets.

CONCLUSION
In this study, we used an ANN to predict the VF using features extracted from 120 s HRV signals, the QRS complex signed area, and the R-peak amplitude 30 s before VF occurrence. The datasets were collected from the popular physiological archive PhysioNet. Although the datasets utilized in this study were relatively small, the performance of the ANN was better using QRS shape features than that of traditional HRV features. This was consistently observed in all machine learning algorithms implemented in this study, which demonstrates the feasibility of using QRS shape features to accurately predict VF onset. This work requires further investigation using a greater number of datasets to confirm the clinical feasibility of our proposed approach. Finally, the results of this study could be used to predict when an ICD will begin ventricular defibrillation.

DATA AVAILABILITY
Publicly available datasets were analyzed in this study. This data can be found here: https://www.physionet.org/.

AUTHOR CONTRIBUTIONS
This manuscript is the intellectual product of the entire team. GT designed the study, wrote the simulation source code and the manuscript, performed the data analysis, and interpreted the results. KL, ES, and H-JH reviewed and revised the whole manuscript based on the simulation results. All authors read and approved the final manuscript.