Detecting and interpreting myocardial infarction using fully convolutional neural networks

As discussed in the previous section, time series classification in a realistic setting has to be able to cope with time series that are so large that they cannot be used as input to a single neural network or that cannot be downsampled to reach this state without losing large amounts of information. At this point two different procedures are conceivable: either one uses attentional models that allow to focus on regions of interest, see e.g. Karim et al (2018a) and (2018b), or one extracts random subsequences from the original time series. For reasons of simplicity and with real-time on-site analysis in mind we explore only the latter possibility, which is only applicable for signals that exhibit a certain degree of periodicity. The assumption underlying this approach is that the characteristics leading to a certain classification are present in every random subsequence. We stress at this point that this procedure does not rely on the identification of beginning and endpoints of certain patterns in the window (Hu et al 2013). This approach can be justified a posteriori with the reasonable accuracies and specificities it achieves. Furthermore, from a medical point of view it is reasonable to assume that ECG characteristics do not change drastically within the time frame of any single recording.

The procedure leaves two hyperparameters: the choice of the window size and an optional downsampling rate to reduce the temporal input dimension for the neural network. As the dataset is not large enough for extensive hyperparameter optimizations we decided to work with a fixed window size of 4 s downsampled to an input size of 192 pixels for each sequence. The window size is sufficiently large to capture at least three heartbeats at normal heart rates.

As discussed in section 3, if we consider a binary classification problem we are dealing with an imbalanced dataset with 80 healthy records in comparison to 127 records diagnosed with MI. Several approaches have been discussed in the literature to best deal with imbalance (Branco et al 2015, Buda et al 2017). Here we follow the general recommendations and oversample the minority class of healthy patients by 2:1.

We refrain from using accuracy as target metric as it depends on the ratio of healthy and infarction ECGs under consideration. As sensitivity and specificity are the most common metrics in the medical context, we choose Youden's J-statistic as target metric for model selection which is determined by the sum of both quantities, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J&={{\rm sensitivity}} + {{\rm specificity}} - 1\nonumber \\ &=\frac{TP}{TP+FN}+\frac{TN}{TN+FP}-1, \label{eq1} \nonumber \end{align} \tag{ 1 }$$

where $TP/FN/FP$ denote true positive/false negative/false positive classification results. Other frequently considered observables in this context include $F_1$ or $F_2$ scores that are defined as combinations of positive predictive value (precision) and sensitivity (recall).

Finally, to obtain the best possible estimate of the test set sensitivity and specificity using the given data, we perform 10-fold cross-validation on the dataset. Its size is comparably small and there are still considerable fluctuations of the final result statistics, even considering the data augmentation via random window selection. These result statistics do not necessarily reflect the variance of the estimator under consideration when applied to unseen data (Isaksson et al 2008) and it is not possible to infer variance information from cross-validation scores by simple means (Bengio and Grandvalet 2004). The given dataset is not large enough to allow a train-validation-test split with reasonable respective sample sizes. Following Shaikhina and Khovanova (2017), we circumvent this problem by reporting ensemble scores corresponding to models with different random initializations without performing any form of hyperparameter tuning or model selection using test set data. Compared to single initializations the ensemble score gives a more reliable estimate of the model's generalization performance on unseen data. For calculating the ensemble score we combine five identical models and report the ensemble score formed by averaging the predicted scores after the softmax layer (Ju et al 2017).

We investigate both convolutional neural networks as well as recurrent neural network architectures. While recurrent neural networks seem to be the most obvious choice for time series data, see e.g. Hüsken and Stagge (2003), convolutional architectures have been applied for similar tasks in early days, see e.g. Waibel et al (1989) for applications in phoneme recognition.

We study different variants of convolutional neural networks inspired by several successful architectures applied in the image domain such as fully convolutional networks (Long et al 2015) and resnets (He et al 2016a, 2016b, Zagoruyko and Komodakis 2016); see appendix A.1 for details. In addition to architectures that are applied directly to the (downsampled) time series data, we also investigate the effect of incorporating frequency-domain input data obtained by applying a Fourier transform to the original time-domain data. We stress again that our approach operates directly on the (downsampled) input data without any preprocessing steps.

For comparison we also consider recurrent neural networks, namely LSTM (Hochreiter and Schmidhuber 1997) cells. We investigate two variants: in the first approach we feed the last LSTM output into a fully connected layer. In the second case we additionally apply a time-distributed dense layer, i.e. a layer with shared weights across all time steps to train the network in addition on a time-series classification task where we adjusted both loss functions to reach similar values. Similar to Hüsken and Stagge (2003) we investigate in this way if the time series predication task improves the classification accuracy.

In table 3 we compare the architectures described in the previous section based on 12-lead data. The comparison is based on cross-validated J-statistics without implying statistical significance of our findings. In the light of very small number of 20 or even fewer patients in the respective test sets, we do not report confidence intervals or similar variance measures as these would be mainly driven by the fluctuations due to the small size of the test set, see also Bengio and Grandvalet (2004), Isaksson et al (2008), Vanwinckelen and Blockeel (2012) and the related discussion in section 4.1.

Table 3. Classification results for different network architectures on 12-lead data.

Model	J-Stat	Sensitivity (Recall)	Specificity	Precision (Positive predicted value)
Fully convolutional	0.827	0.933	0.897	0.936
Resnet	0.828	0.925	0.903	0.940
Fully convolutional(freq)	0.763	0.902	0.860	0.913
Resnet(freq)	0.656	0.870	0.786	0.869
LSTM mode (final output)	0.743	0.910	0.833	0.899
LSTM (final output $+$ pred.)	0.742	0.914	0.828	0.897

The fully convolutional architecture and the resnet achieve similar performance applied to time-domain data. In contradistinction to an earlier investigation (Wang et al 2017) that favored the fully convolutional architecture, a ranking of the two convolutional architectures is not possible on the given data. Interestingly, the convolutional architectures perform better applied to raw time-domain data than applied to frequency-domain data. In this context it might be instructive to investigate also other transformations of the input data such as wavelet transformations as considered in Sharma et al (2015) and Sharma and Sunkaria (2018).

Both convolutional architectures show a better score than recurrent architectures. This can probably be attributed to the fact that we report just results with standard LSTMs and do not investigate more advanced mechanisms such as most notably an attention mechanism, see e.g. Song et al (2017), Karim et al (2018a), (2018b) for recent developments in this direction. Training the recurrent neural network jointly on a classification task as well as on a time series prediction task, see also the description in appendix A.2, did not lead to an improved score, whereas a significant increase was reported by Hüsken and Stagge (2003), which might be related to the small size of the dataset.

In the following sections we analyze particular aspects of the classification results in more detail. All subsequent investigations are carried out using the fully convolutional architecture, which achieved the same performance as the best performing resnet architecture with a comparably much simpler architecture. If not noted otherwise we use the default setup of 12-lead data.

5.2.1. MI localization and training procedure

As described in section 3 we distinguish the aggregated subdiagnosis classes aMI and iMI. Here we examine the classification performance of a model that distinguishes these subclasses rather than training just on a common superclass MI. We can investigate a number of different combination of either training/evaluating with or without subdiagnoses as shown in table 4.

Table 4. Classification performance on subdiagnoses (with fully convolutional architecture) on 12-lead data. Train MI (train aMI$+$iMI) refers to training disregarding (incorporating) subdiagnoses and train aMI/iMI to training using just a particular subdiagnosis; analogously for eval.

Data	J-Stat	Sensitivity (Recall)	Specificity	Precision (Positive predicted value)
Cardiologists aMI (Willems et al 1991)	0.857	0.874	0.983	—
Cardiologists iMI (Willems et al 1991)	0.738	0.749	0.989	—
Train MI eval MI	0.827	0.933	0.897	0.936
Train MI eval aMI	0.877	0.980	0.897	0.884
Train MI eval iMI	0.789	0.894	0.896	0.879
Train aMI eval aMI	0.880	0.919	0.961	0.950
Train iMI eval iMI	0.689	0.810	0.879	0.849
Train aMI$+$iMI eval MI	0.788	0.912	0.876	0.947
Train aMI$+$iMI eval aMI	0.846	0.966	0.881	0.906
Train aMI$+$iMI eval iMI	0.741	0.861	0.879	0.897

Both for models trained on unspecific infarction and for models trained using subdiagnosis labels, the performance on the iMI classification task turns out to be worse than the score achieved for aMI. The most probable reason for this is that anterior myocardial infarctions show typical signs in most of the Wilson leads because of the proximity of the anterior myocard to the anterior chest wall. For the more difficult task of iMI classification, the model seems to profit from general MI data during training, as a model trained on generic MI achieves a higher score on aMI classification than a model trained specifically on aMI classification only. The converse is true for the simpler task of aMI classification.

Interestingly, the model trained without subdiagnoses reaches a slightly higher score both for unspecific MI classification as well as for classification on subdiagnoses aMI/iMI only, which might just be an effect of an insufficient amount of training data. In any case, we restricted the rest of our investigations on the model trained disregarding subdiagnoses.

In figure 1 we show the confusion matrix for the model that is trained and evaluated on the subdiagnoses aMI and iMI in addition to healthy control. The confusion matrix underlines the fact that the model is able to discriminate between the aggregated subdiagnoses whose assignments were motivated by medical arguments, see section 3, and represents an a posteriori justification for this choice.

Zoom In Zoom Out Reset image size

The reported score for models evaluated on subdiagnoses allows a comparison to human performance on this task. We base this comparison on a study (Willems et al 1991) that assessed the human classification performance for different diagnosis classes based on a panel of eight cardiologists. Here we only report the combined result and refer to the original study for individual results. The most appropriate comparison is the model trained on general MI and evaluated on subdiagnoses aMI and iMI. However, it turns out that irrespective of the training procedure the algorithm achieves a slightly higher score on aMI classification and a considerably higher score on iMI classification and we see it therefore justified to claim at least human-level performance on the given classification task. We restrain from drawing further conclusions from this comparison as it depends on the precise performance metric under consideration, the fact that not the same datasets were used in both studies and differences in subdiagnosis assignments. The claim of superhuman performance on this task would certainly require more thorough investigations in the future.

5.2.2. Comparison to literature approaches

Our data and channel selection strategy, see section 3, was carefully chosen to reflect the requirements of a clinical application as closely as possible. In addition, considering the comparably small size of the dataset, a careful cross-validation strategy is of utmost importance, see section 4.1. Unfortunately, most literature results do not report cross-validated scores or introduce data leakage in their cross-validation procedures. This can happen for example by sampling from beat-level segmented signals or most commonly by sampling based on ECGs rather than patients; see Sharma and Sunkaria (2018) for a detailed discussion. Both cases lead to unrealistically good performance estimates as the classification algorithm can in some form adapt to structures in the same ECG or structures in a different ECG from the same patient during the training phase.

We refrain from presenting results for the latter setups as they do not allow to disentangle a model's classification performance from its ability to reproduce already known patterns. Therefore, we only include a comparison to the most recent works Reasat and Shahnaz (2017) and Sharma and Sunkaria (2018) that are to our knowledge the only works where a cross-validated score with sampling on patient level, in the literature also termed subject-oriented approach (Sharma and Sunkaria 2018), is reported. To ensure comparability with Reasat and Shahnaz (2017) and Sharma and Sunkaria (2018) we replicate their setup as closely as possible and modify our data selection to include limb leads and in addition to healthy records only all genuine iMI ECGs. In this case our approach shows not only superior performance compared to literature results, see table 5, but does unlike their algorithms operate directly on the input data and does not require preprocessing with appropriate input filters.

Table 5. Comparison to literature results (with fully convolutional architecture). Results from this work marked by ^a.

Benchmark	J-Stat	Sensitivity (Recall)	Specificity	Precision (Positive predicted value)
Wavelet transform $+$ SVM (Sharma and Sunkaria 2018)	0.583	0.790	0.793	0.803
CNN (Reasat and Shahnaz 2017)	0.694	0.853	0.841	—
Limb leads $+$ inferior MIa	0.773	0.874	0.900	0.932

We replicated the above setup to demonstrate the competitiveness of our approach, but for a number of reasons we are convinced that the scores presented in tables 3 and 4 are the more suitable benchmark results: firstly, from a clinical point of view 12-lead ECGs are the default choice and the algorithm should be fed with the full set of eight non-redundant channels. Secondly, the restriction to include only the first infarction ECG per patient is arguably more suited for the application of the clinically most relevant problem of classifying acute myocardial infarctions, see the discussion in section 3. Finally, from a machine learning perspective it is beneficial to include all subdiagnoses for training allowing to adapt to general patterns in infarction ECGs and to only evaluate the trained classifier on a particular subdiagnosis of interest. For the case of aMI this procedure leads to an improved score, see the discussion of table 4.

5.2.3. Channel selection

A general challenge remains the topic of interpretability of machine learning algorithms and in particular deep learning approaches that is especially important for applications in medicine (Cabitza et al 2017). In the area of deep learning, there has been a lot of progress in this direction (Simonyan et al 2013, Sundararajan et al 2017, Montavon et al 2018). So far most applications covered computer vision whereas time series data in particular did only receive scarce attention. Interpretability methods have been applied to time series data in Wang et al (2017) and ECG data in particular in Teijeiro et al (2018). A different approach towards interpretability in time series was put forward in Siddiqui et al (2018).

As an exploratory study for the application of interpretability methods to time series data we investigate the application of attribution methods to the trained classification model. This allows investigating on a qualitative level if the machine learning algorithm uses similar features as human cardiologists. Our implementation makes use of the DeepExplain framework put forward in Ancona et al (2018). For neural networks with only ReLU activation functions it can be shown (Ancona et al 2018) that attention maps from '${{\rm gradient}}\times {{\rm input}}$' (Shrikumar et al 2016) coincide with attributions obtained via the $\newcommand{\e}{{\rm e}} \epsilon$-rule in LRP (Bach et al 2015). Even though we are using ELU activation functions the attribution maps show only minor quantitative differences. The same applies to the comparison to integrated gradients (Sundararajan et al 2017). For definiteness, we focus our discussion on '${{\rm gradient}}\times {{\rm input}}$'. Different from computer vision, where conventionally attributions of all three color channels are summed up, we keep different attributions for every channel to be able to focus on channel-specific effects. We use a common normalization of all channels to be able to compare attributions across channels.

We stress that attributions are inherent properties of the underlying models and can therefore differ already for models with different random initializations in an otherwise identical setup. If we aim to use it to identify typical indicators for a classification decision as a guide for clinicians a more elaborate study is required. For simplicity, in this exploratory study we focus on a single model rather than the model ensemble. By visual inspection we identified the most typical attribution pattern for MI among examples in the batch that occurred shortly after the infarction. Prototypical outcomes of this analysis are presented below.

Figure 2 shows examples of the interpretability analysis of selected channels of two MI ECGs that were correctly classified. As in the clinical context ECGs are always considered in the context of the full set of 12 channels (if available), the complete set of channels is shown in figure 3 in the appendix. Take note that the attributions show a high consistency over beats of one ECG, even if a significant baseline shift is present. There is also a reasonable consistency with regard to similar ECG features exhibited by other patients which are not shown here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

ECG A is taken from a 74 year-old male patient one day after the infarction took place. A coronary angiography performed later confirmed an aMI. ECG B is taken from a 68 year-old male patient one day after the infarction which was later confirmed to be in inferior localization. ECGs A and B are listed as s0021are and s0225lre in the PTB dataset.

Signs for ischemia and infarction are numerous and of variable specificity (Tewelde et al 2017). Highlighted areas coincide with established ECG signs of MI. These are typically found between and including the QRS complex and the T wave, as this is when the contraction and consecutive repolarization of the ventricles take place. ST-segment-elevation (STE) is the most important finding in MI ECGs and diagnostic criterion for ST-elevation myocardial infarction (STEMI) (Thygesen et al 2012). This sign (though not formally significant in every case) and corresponding high attributions can be found in both example ECGs at the positions marked ${\bf a}$, ${\bf e}$, and ${\bf g}$. At the same time in another channel (position ${\bf c}$) there is no STE and the attribution is consequently inverted. The attribution at position ${\bf a}$ also coincides with pathological Q waves, which also occur in some infarction ECGs. T wave inversion, another common sign for infarction, can be found at position ${\bf d}$. Some other attributions of the model are less conclusive. Although attributions at positions ${\bf b}$ and ${\bf f}$ fall in the T and/or U waves, that is regions that are relevant for detection of infarction, it is unclear why they influenced the decision against infarction.

Note that the highlighted areas do not necessarily align perfectly with what clinicians would identify as important. For example, for a convolutional neural network to detect an ST-elevation, it must use and compare information from before and after the QRS complex, which most likely results in high attributions to the QRS complex itself and its immediate surrounding rather than to the elevated ST-segment.

Comparing the overall visual impression of the attributions across all channels (see figure 3), the model seems to attribute more importance to the Wilson leads in ECG A (anterior infarction) and more importance to the limb leads in ECG B (inferior infarction). This is also where clinicians would expect to find signs of infarction in these cases.

Attributions are inherently model-dependent and as a matter of fact the corresponding attributions show quantitative and in some cases even qualitative differences. However, across different folds and different random initializations the attribution corresponding to the STE was always correctly and prominently identified. This is a very encouraging sign for future classification studies on ECG data based on convolutional methods, in particular in combination with attribution methods. A future study could put the qualitative finding presented in this section on a quantitative basis. This would require a segmentation of the data, possibly using another model trained on an annotated dataset as no annotations are available for the PTB dataset, and statistically evaluating attribution scores in conjunction with this information. In this context it would be interesting to see if different classification patterns arise across different models or if they can at least be enforced as in Ross et al (2017).

We generally use ELU (Clevert et al 2015) as activation function both for convolutional as well as fully connected layers without using batch normalization (Ioffe and Szegedy 2015), which was reported to lead to a slight performance increase compared to the standard ReLU activation with batch normalization (Mishkin et al 2017). In the architectures with fully connected layers we apply dropout (Srivastava et al 2014) at a rate of 0.5 to improve the generalization capability of the model. We initialize weights according to He et al (2015). Note that in contrast to the case of 2D data a max pooling operation only reduces the number of couplings by a factor of 2 rather than 4, which is then fully compensated by the conventional increase of filter dimensions by 2 in the next convolutional layer. To achieve a gradual reduction of couplings we therefore keep the number of filters constant across convolutional layers. We study the following convolutional architectures that are also depicted in figure A1.

• (1)  
  A fully convolutional architecture (Long et al 2015) with a final global average pooling layer.
• (2)  
  A resnet-inspired (He et al 2016a, 2016b, Zagoruyko and Komodakis 2016) architecture with skip-connections.

Zoom In Zoom Out Reset image size

We investigate the impact of including frequency information obtained via a fast Fourier transformation with $N_{\rm FFT}=2^{\lceil{\log_2 (d)}\rceil}$, where $d=192$ denotes the sequence length after rescaling. The $N_{\rm FFT}/2+1$ independent components are used as frequency-domain input data with otherwise unchanged network architectures.

As alternative architecture we investigate recurrent neural networks, namely LSTM (Hochreiter and Schmidhuber 1997) cells. We investigated stacked LSTM architectures but found no significant gain in performance. However, even for a single RNN cell, in our case with 256 hidden units, different training methods are feasible.

• (1)  
  In the first variant we feed the last LSTM output into a fully connected softmax layer.
• (2)  
  In the second variant we additionally apply a time-distributed fully connected layer, i.e. a fully connected layer with shared weights for every timestep, and train the network to predict the next element in a time series prediction task jointly with the classification task. Here we adjusted both loss functions to reach similar values. Similar to Hüsken and Stagge (2003) we investigate in this way if the time series predication task improves the classification accuracy.

During RNN training we apply gradient clipping.