Automatic detection of T wave alternans using tensor decompositions in multilead ECG signals

Griet Goovaerts; Bert Vandenberk; Rik Willems; Sabine Van Huffel

doi:10.1088/1361-6579/aa7876

1. Introduction

Sudden cardiac death caused by cardiac arrhythmia is one of the major causes of death worldwide (Wellens et al 2014). Since many arrhythmias can be reversed with an electrical shock from an (implantable) cardiac defibrillator (ICD), it is important to identify patients at risk for such life-threatening events in order to identify candidates for a preventive implantation of an ICD (Priori et al 2015). T wave alternans (TWA) can be derived from electrocardiogram (ECG) and is a promising feature for risk stratification. Studies have shown that the presence of TWA is correlated with the risk of arrhythmia in patients with ischemic cardiomyopathy after previous myocardial infarctions and long QT syndrome (Shimizu and Antzelevitch 1999, Ikeda et al 2002, Chow 2006). T wave alternans is defined as a periodic beat-to-beat variation in the amplitude of the T wave. Large and small T waves alternate in an ABABAB-pattern. Since the T wave represents the repolarization of the ventricles, it is linked with the repolarization characteristics of the heart. When the variation between subsequent T waves is large enough, TWA can be detected by visually inspecting the ECG signal. In cases where the difference is only a few microvolts or when dealing with long ECG signals this is however not feasible. Therefore (semi-)automatic methods are used that can reliably detect TWA in the ECG signal.

Martínez and Olmos (2005) give an overview of commonly used algorithms for (semi-)automatic T wave alternans detection. Most methods start with a preprocessing stage followed by the actual TWA detection. Examples of commonly used methods are the correlation method (Burattini et al 1999), modified moving average method (Nearing and Verrier 2002) or spectral method (Smith et al 1988). Most methods decompose the ECG signal to a beat-to-beat time series that contains the T wave characteristics. The actual detection of T wave alternans is then done on this decomposed signal. The majority of TWA detection methods are designed for use with single channel ECG signals. Multichannel ECG signals have however several advantages over single lead ECGs. When the ECG is measured over different channels it contains information about the behavior of the heart in multiple spatial dimensions, which can prompt a more comprehensive understanding about certain cardiac disorders. Furthermore when there is noise present in one or more channels, the use of a multilead ECG signal allows (partial) reconstruction of the original signal. This is obviously impossible if there is only one lead present. Tensors are a natural way of dealing with multichannel signals. They have been used widely in chemometrics and psychometrics since they have many properties that are superior to standard matrix methods (Kolda and Bader 2008). Recently there have been applications of tensors in biomedical signal processing, but to date they have mainly been used in the processing of neural signals (Hunyadi et al 2016, Lee et al 2007, Hunyadi et al 2014, Zink et al 2016). Their use in ECG processing is limited to specific applications such as fetal ECG extraction and classification of myocardial infarction (Debals et al 2016, Padhy and Dandapat 2017). We wanted to show that tensors are a robust and reliable way of detecting and quantifying T wave alternans.

The objective of this paper is to present a novel TWA detection method based on tensor decompositions. The method is fully automated and both detects and quantifies the amount of T wave alternans in a multilead ECG signal. The signal is segmented and reorganized in a tensor. By decomposing the tensor and analyzing the resulting factor vectors, TWA can be detected and the amount of TWA can be estimated. Preliminary results were presented in Goovaerts et al (2014) and Goovaerts et al (2015b). The method has now been fully validated on both artificial and clinical signals. Furthermore the use of the more general PARAFAC2 tensor decomposition method makes application in realistic situations possible. In the next section the datasets on which the method has been evaluated are described. Section 3 describes the different steps necessary to detect T wave alternans, including a description of the tensor decomposition methods. The results are summarized in section 4, followed by discussion and conclusion.

2. Data

The method has been evaluated on three different datasets: a simulation dataset, the TWA challenge database available on Physionet and a dataset from the University Hospitals Leuven.

2.1. Simulation data

Artificial ECG and noise signals were generated using the OSET toolbox (Sameni 2012). The artificial signals consist of the 12 standard ECG channels and are 30 s long. The sampling frequency is 250 Hz. The signals are simulated by modeling a cardiac dipole and linearly projecting it onto the body surface using the Dower transform (Dower et al 1988). T wave alternans was simulated by introducing a constant offset to each odd T wave as demonstrated in Clifford et al (2008). The level of TWA was varied between 5 μV and 50 μV. After all signals were generated, the differences between subsequent T waves were manually checked to ensure the TWA level was consistent with the simulated levels. Fifty different patients are simulated by calculating 50 Dower transformation matrices based on signals of the PTB database available on Physionet (Bousseljot et al 1995, Goldberger et al 2000). In order to work with realistic ECG signals and to be able to test the robustness to noise, realistic ECG noise was added to the signals using the noise generator provided with (Sameni 2012). The noise consists of a mixture of equal amounts of baseline wander, muscle artifacts and electrode movements. A separate noise signal was generated for each synthetic signal and superposed onto the signal. The SNR value of the signals was varied to construct two sets with respectively moderate noise ( ${\rm SNR} = 15~{\rm dB}$ ) and high noise ( ${\rm SNR} = 5~{\rm dB}$ ). Figure 1 illustrates the noise level by showing a signal from one patient in all three conditions (without noise, moderate noise, heavy noise).

**Figure 1.** Example of a simulated signal from the same patient without noise (left), moderate noise (middle) and heavy noise (right).
Download figure:
Standard image High-resolution image

2.2. Physionet dataset

The TWA database available on Physionet (Goldberger et al 2000) was assembled for the computing in cardiology challenge of 2008 (Moody 2008). It contains 100 ECG signals sampled at 500 Hz. Most signals are 12-lead ECG channels, some signals only contain two or three channels. The duration of the signals is (approximately) two minutes. The records were collected from different databases, and the final dataset contains both real-life and artificial ECG signals. The signals are ranked according to the level of T wave alternans present in the ECG.

2.3. University Hospitals Leuven dataset

The dataset was obtained in the University Hospitals Leuven (Belgium), and includes 33 Holter ECG signals from 9 subjects. Selection of subjects was done based on the results of a clinical TWA test (Cambridge Analytic Spectral Method). Only subjects with multiple consistent positive or negative tests were included in the dataset. During the test, subjects were asked to cycle on a stationary bike. By increasing the workload, their heart rate is increased from resting rate to 110 beats per minute. The dataset included 13 signals from 4 subjects with positive clinical TWA tests (referred to as TWA group) and 20 signals from 5 subjects with negative clinical results (referred to as control group). All ECG signals are recorded with a sampling frequency of 256 Hz. Most records have three ECG channels (30/33 records), some only include two channels (3/33 records). The signals are all between 14 and 23 min long. The study was approved by the Ethical Committee of the University Hospitals Leuven.

3. Method

The complete TWA detection method consists of five steps. After preprocessing and T wave segmentation, a tensor containing all T waves is constructed. This tensor is then decomposed and the decomposition result is analyzed to detect and quantify T wave alternans.

3.1. Preprocessing

The preprocessing stage is necessary to remove noise from the signal so that it does not affect the results of the TWA detection. When applying noise removal methods, it is essential that the shape and amplitude of the T wave are not altered since this could introduce errors in the final result. We consider baseline drift and high frequency noise from for example muscle artifacts as major noise sources in ECG signals. To remove baseline wander a method that removes the baseline by reducing quadratic variation is used (Fasano and Villani 2014). After that, the signal is filtered with a low-pass filter with a cut-off frequency of 60 Hz to remove high frequency noise. Since the T wave is typically contained in a frequency band ranging from 2 to 35 Hz (Burattini and Giuliani 2013), this should not induce significant changes in morphology of the T wave on a macroscopic scale. Finally the signals are normalized to facilitate comparison between different signals.

3.2. T wave segmentation

In order to be able to detect transient T wave alternans (i.e. TWA that is only present in certain segments of the ECG), a moving window of 50 heartbeats is used to analyze the signals. T wave segmentation is preferred over T wave delineation since the latter is sensitive to noise. T waves are segmented by selecting a window with length wl after each detected R peak. Since the length of the QT interval will change with a changing heart rate, the start point st and length of the windows wl are dependent on the mean RR interval length over the current analysis window of 50 heartbeats ( $\overline{RR}$ ). For st we make a distinction between three possible $\overline{RR}$ ranges:

(i)
$\overline{RR} \leqslant 0.6~s$ : st = R peak + 100 ms
(ii)
$0.6~s <\overline{RR} <1.1~s$ : st = R peak + 150 ms
(iii)
$\overline{RR} \geqslant 1.1s$ : st = R peak + 250 ms

The window length wl is $0.4\sqrt{\overline{RR}}$ . The interval boundaries have been determined experimentally to ensure accurate T wave segmentation and largely correspond to other values found in literature (Burattini et al 1999).

3.3. Tensor construction

Tensors are by definition multidimensional arrays with a number of dimensions greater than two. An ECG signal naturally is a matrix with dimensions channels × time. To transform this 2D matrix to a multidimensional tensor we can either add extra data (such as the frequency content of the signal) or restructure the raw signal in a 3D manner. Since the objective is T wave alternans detection, the tensor should be constructed in such way that the changes in T waves over different heartbeats are maximally emphasized. Therefore we will construct a T wave tensor: a three-dimensional tensor that consists of the T waves of all channels and heartbeats. The T waves of heartbeat one will form the first slice of the tensor, the T waves of the second heartbeat form the second slice, etc. This way, one tensor is constructed for each ECG signal. Each tensor has three dimensions, channels × time × heartbeats. If the ECG signal contains N channels and M heartbeats, the resulting tensor will have size ( $N \times wl \times M$ ). Since we use a sliding window of 50 heartbeats to analyze the signals, M is here fixed to 50. Note that both the second and third dimensions contain temporal information: the second dimension contains the signal of individual T waves while the third dimension shows a more long-term view of the different heartbeats. Figure 2 illustrates the construction of the tensor.

**Figure 2.** Construction of the T wave tensor: the T waves of all channels are stacked beat-by-beat in a 3D manner. The result is a tensor with dimensions *channels* × *time* × *heartbeats*.
Download figure:
Standard image High-resolution image

3.4. Tensor decomposition

3.4.1. Canonical polyadic decomposition (CPD).

Many types of tensor decomposition methods exist which will break down the tensor in to different factors. Canonical polyadic decomposition (CPD) decomposes a multidimensional tensor X in a (minimal) sum of R rank-1 terms (Kolda and Bader 2008). One of the main advantages of CPD over standard matrix decomposition methods is that it is unique up to permutation and scaling under mild conditions. If X has three dimensions CPD can be expressed as:

$\begin{align} \displaystyle X = \sum_{r=1}^R A_r \circ B_r \circ C_r + \epsilon. \nonumber \end{align} \tag{ 1 }$

The number of rank-1 terms R is the rank of the decomposition, $A_r, B_r$ and C_r are the factor vectors for each dimension and epsilon the model error. Since in this application the interest lies in the main variation between subsequent T waves, only 1 component is extracted and therefore R is restricted to 1. A_r, B_r and C_r will then be 3 vectors of size $N \times 1$ , $wl \times 1$ and $M \times 1$ respectively that correspond to the three dimensions of the original tensor, channels, time and heartbeats.

CPD results in one loading vector for each tensor dimension. If the tensor is low-rank, the CPD model is accurate and these 3 vectors will capture the main variations present in the tensor. When the ECG signal is relatively noisefree and the T waves are perfectly aligned in the third dimension, this assumption is met and CPD will lead to correct results. When there is however a variation in, for example, the timing of subsequent T waves (as can be the case with a changing heart rate), the tensor cannot be decomposed in rank-1 components and CPD will lead to inaccurate results. In those cases it is better to use a more general tensor decomposition method such as PARAFAC2.

3.4.2. PARAFAC2.

PARAFAC2 is a tensor decomposition method that is a variation of the standard CPD (Bro et al 1999). The main difference is that PARAFAC2 is less restrictive than CPD in the sense that it allows variations in the factor vectors of one dimension. This is especially useful when one of the factors contains a time shift. In this application, such time shift arises if the T waves are not exactly aligned in time, for example when the heart rate differs and the QT interval length changes. The difference is shown in figure 3, which shows schematic representations of CPD (figure 3(a)) and PARAFAC2 (figure 3(b)). If X_k is one frontal slice of a tensor with the T waves of one heartbeat in all channels, this slice can be modeled as:

$\begin{align} \displaystyle X_m = \sum_{r=1}^R A_r \circ B_{rm} \circ C_r + \epsilon. \label{eq:PARAFAC2} \nonumber \end{align} \tag{ 2 }$

**Figure 3.** Illustration of tensor decomposition methods. (a) Canonical polyadic decomposition. (b) PARAFAC2.
Download figure:
Standard image High-resolution image

The subscript m implies that there is an individual loading vector $B_{rm}$ for each frontal slice of the tensor, effectively making B_r a loading matrix with dimensions $wl \times M$ . Each row of this loading matrix corresponds to the temporal profile of a T wave in one heartbeat, taking into account the possible time shifts in T waves due to heart rate changes. M is thus the total number of heartbeats contained in one tensor, and wl the length of the T wave segmentation window. Equation (2) will have a unique solution when the cross-product $B_rm^TB_rm$ remains constant for all m. This will be the case if the different rows of B are only shifted in time. Under normal circumstances, e.g. when the ECG signal does not contain an excessive amount of noise, this assumption is met and PARAFAC2 will accurately decompose the ECG signal.

3.5. Tensor decomposition results

The result of the tensor decomposition is a collection of 3 factor vectors (or 2 factor vectors and 1 matrix in the case of PARAFAC2), which give information about the structure of the signal in the three tensor dimensions. All three components give valuable information about the T wave characteristics in time and space, but only one will be used for the actual detection of T wave alternans. The first mode-1 factor vector corresponds to the spatial dimension. This vector shows the variation of the T waves over the different channels. The mode-2 factor vector (or matrix for PARAFAC2) shows the time course of the T waves. If there is a time shift present this will manifest itself in the different rows of the PARAFAC2 factor matrix. In the case of CPD this vector simply shows the average T wave in the signal. The mode-3 factor vector expresses the differences in the T wave in subsequent heartbeats. Only this vector is used in the next step for further analysis.

3.6. Detection of T wave alternans

The actual detection of T wave alternans is done using the third factor vector as explained in the previous paragraph. If a signal contains TWA, the typical ABABAB-pattern that is present in the T wave amplitude will be captured by the tensor decomposition as an alternating time series in the third dimension. To detect TWA the number of consecutive turning points of C_r, the loading vector of the heartbeats dimension, are first detected. Turning points are defined here as local minima or maxima, e.g. points where the derivative of the factor vector changes sign. Since noise also causes (random) variations in the T wave amplitude, a threshold on the minimum number of consecutive turning points is introduced to avoid false positive detections. The threshold is set to 10, since a preliminary study showed that this value leads to good results with real-life signals (Goovaerts et al 2015b). The TWA magnitude is then determined by summing the differences between all turning points.

Let TP be the vector of turning point indices (taking into account the minimal number of 10 consecutive turning points) with length L, and C_r the third loading vector in the result of the tensor decomposition. The total TWA magnitude in a window can then be calculated as:

$\begin{align} \displaystyle \label{eq:TWA} {\rm TWA} = \sum_{i=2}^{L} \vert C_r(TP(i)) - C_r(TP(i-1))\vert. \nonumber \end{align} \tag{ 3 }$

4. Results

This section presents the results of TWA detection on different types of data. In all cases a comparison between CPD and PARAFAC2 is done in order to assess the differences between both methods. First the detailed results for three case studies (two signals containing TWA, one with and one without T wave shift, and one signal without TWA) are qualitatively described. Next the results for the three datasets considered are discussed. All experiments were performed using MATLAB R2015b. Tensorlab, a MATLAB-based toolbox containing many numerically optimized methods for tensor calculations and structured data fusion, is used for the tensor construction and CPD decomposition (Vervliet et al 2016). PARAFAC2 is implemented in the N-way toolbox (Andersson and Bro 2000).

4.1. Case study

We examined three sample signals taken from the Physionet database. The first signal (twa11) did not present TWA, the second signal (twa99) contained a significant amount of T wave alternans and the third signal (twa01) contained both TWA and shifted T waves caused by changes in heart rate. This way we can evaluate both the effect of the presence of TWA and a shift in T wave position for both methods. All signals consist of 12 channels and have been cropped to 50 heartbeats. Figure 4 shows a 10 s excerpt of the signals. None of the signals contains obvious noise, meaning the results are not influenced by signal quality. The complete signals (30 min, 12 channels) can be found on https://physionet.org/pn3/twadb/. After tensorization and tensor decomposition of all three signals (using both CPD and PARAFAC2), the resulting loading vectors could be examined. The results are shown in figure 5.

**Figure 4.** Examples 10 s of all three signals used in the case study, taken from the Physionet TWA database. (a) Sample of the signal without TWA (*twa11*). (b) Sample of the signal with TWA (*twa99*). (c) Sample of the signal with TWA and T wave shift (*twa01*).
Download figure:
Standard image High-resolution image

**Figure 5.** Values of factor vectors A, B and C for signal without TWA (top row), a signal with TWA (middle row) and a signal with TWA and shifted T waves (bottom row). The results for CPD are plotted in red, the results for PARAFAC2 in black. The three different columns show respectively the spatial, temporal and heartbeats feature vectors. (a) Spatial vector (no TWA). (b) Temporal vector (no TWA). (c) Heartbeats vector (no TWA). (d) Spatial vector (TWA). (e) Temporal vector (TWA). (f) Heartbeats vector (TWA). (g) Spatial vector (TWA + shift). (h) Temporal vector (TWA + shift). (i) Heartbeats vector (TWA + shift).
Download figure:
Standard image High-resolution image

Figures 5(a), (d) and (g) show the spatial vector for the signal without TWA, with TWA and with TWA and T wave shift respectively which show the sign and magnitude of the T wave in different channels. From figure 5(a) we can for example conclude that T waves in channel 9 will have a larger amplitude than those in channel 12. The temporal vectors, shown in figures 5(b), (e) and (h) correspond to the shape of the T wave in time. Here the difference between CPD and PARAFAC2 is obvious, especially in figure 5(h). The result for PARAFAC2 will be a matrix, and each row of the matrix (corresponding to the temporal vector for each heartbeat) is plotted on top of each other. The presence of a large T wave shift in the third signal is evident. While CPD succeeds in capturing the mean shape of the T wave, the difference in the timing of the subsequent waves is discarded. This information is contained in the feature matrix obtained by PARAFAC2. The other signals do not contain a large shift, as is demonstrated by the small variation in the PARAFAC2 results of figures 5(b) and (e). Hence, the results of CPD and PARAFAC2 will be much more similar. T wave alternans can be detected by analyzing the third loading vector. Figures 5(c), (f) and (i) show this loading vector for all analyzed signals. On figure 5(c) we see that, although there is a certain variation in the T wave magnitude over different heartbeats, no clear pattern is visible. The typical ABABAB-pattern is however clearly visible on figures 5(f) and (i). Note that in figure 5(f), there is only T wave alternans present in the first part of the signal. The total TWA magnitude can then be calculated as defined in equation (3). Note that the results for PARAFAC2 and CPD are almost identical in figures 5(c) and (f), while the results in figure 5(i) differ significantly. This is caused by the large T wave shift (see figure 5(i)).

4.2. Artificial signals

We first analyzed the results of both proposed methods on a clean signal with varying levels of T wave alternans. Then we investigated the influence of a temporal shift in the T waves. Finally the robustness of the methods was tested by changing the noise levels. All results are summarized in figure 6. All results shown are averaged over the 50 simulated patients, with the variance among all patients indicated with the shaded areas.

**Figure 6.** Results for CPD (red) and PARAFAC2 (black) for four types of artificial signals: clean signals with varying amount of TWA (a), clean signals ( ${\rm TWA} = 25$ μV) with changing T wave shift (b), artificial signals with a moderate (c) and high (d) noise level. The circles represent the mean values over all 50 simulated patients, the shaded areas the variations in the results. (a) Simulated vs estimated TWA levels for artificial signals without noise. (b) Estimated TWA levels for different values of T wave shift. (c) Simulated vs estimated TWA levels for artificial signals with moderate noise. (d) Simulated vs estimated TWA levels for artificial signals with high noise.
Download figure:
Standard image High-resolution image

The results for a signal without noise are presented in figure 6(a). The simulated TWA level is varied between 5 μV and 50 μV. There is a clear correlation between the simulated and estimated TWA magnitude, both for CPD and PARAFAC2. There is very little variance between results for all patients: the mean differences in estimated TWA levels over all patients is 3 μV and does not change significantly with the simulated TWA level. An artificial T wave shift is then introduced by moving the T wave segmentation window a random number of samples (figure 6(b)). The maximal shift is varied between 0 and 50 ms and the TWA level is fixed to 25 μV. The results are presented as the ratio between the estimated TWA and simulated TWA, which should be 1 for a good detection. The results for PARAFAC2 stay approximately constant (between 0.98 and 1.01), even when the T wave shift is maximal. CPD does not handle the increasing T wave shift well: for a very small shift (between 0 and 10 ms) the result is correct, but for larger shifts the results deteriorate quickly. For a maximal T wave shift of 50 ms the average CPD result seems good, but the large variance among different patients indicate that the obtained results vary greatly and are thus not reliable. When introducing noise, the correlation between the simulated and estimated TWA levels is maintained, although the estimated CPD levels are distinctly smaller figure (6(c)). For a simulated TWA level of 50 μV, the average result obtained with CPD is 39 μV, while the result of PARAFAC2 is more accurate (43 μV). The results show a larger spread for lower TWA levels, which could be expected since the noise is more likely to mask the TWA in these cases. Both methods show less precise results under high noise level (see figure 6(d)): while the estimated TWA magnitudes on average increase with increasing simulated TWA, there are clearly more variations when compared with figures 6(a) and (c). Here PARAFAC2 seems to underperform compared to CPD, especially for low TWA levels.

4.3. Physionet database

To our knowledge, the Physionet database is the only publicly accessible labeled TWA database. This makes it ideal for use as a benchmarking tool to make comparison with existing methods possible. This can be done by ranking the results from the TWA detection in order of magnitude, and comparing this ranking with the reference ranking by calculating the Kendall rank correlation coefficient between the two. A Kendall coefficient that is greater than 0.436 is considered significant ( $p > 0.99$ ). Both an order using CPD and PARAFAC2 have been computed. The Kendall coefficients obtained were both significant, the CPD coefficient being 0.79 and the PARAFAC2 coefficient 0.87. Table 1 compares the tensor-based scores with some methods found in literature and the best result of the original TWA challenge. A coefficient of 0.87 and 0.79, achieved using respectively PARAFAC2 and CPD, would have been the fourth and sixth highest rank obtained during the challenge, and are both higher than the values for some of the most widely used methods found in literature.

Table 1. Kendall coefficient scores obtained by comparing the rankings from different methods found in literature and the two proposed tensor-based methods with the reference ranking for the Physionet database.

Method	Kendall coefficient
CPD	0.7851
PARAFAC2	0.87
Periodic component analysis (Monasterio et al 2010)	0.766
Modified moving average method (Nearing and Verrier 2002)	0.73
Spectral method (Smith et al 1988)	0.416
Best challenge (Sieed and Hasan 2008)	0.911

Figure 7 shows a comparison between the reference ranking of the signals and the ranking obtained by respectively CPD (figure 7(a)) and PARAFAC2 (figure 7(b)). There is a clear correlation for both methods, which is reflected in the high Kendall coefficient. The reference ranking and PARAFAC2 ranking correspond largely, and for most signals there is only a small difference between the reference score and the PARAFAC2 score. A remarkable exception is signal twa63, which has an intermediate amount of TWA according to the reference ranking (position 52) but where no TWA is detected using the tensor-based approach (position 3 with both CPD and PARAFAC2). Further inspection of this signal showed that this signal is taken from the Sudden Cardiac Death Holter Database (Greenwald 1986), a database containing signals from patients with sustained ventricular tachyarrhythmia. In this signal in particular, the timing of the T waves is considerably later than in other signals, leading to a suboptimal segmentation and tensor construction. This, combined with the presence of the ventricular tachyarrhythmia, most likely causes the incorrect result in this case.

4.4. Clinical dataset from the University Hospitals Leuven

When analyzing the clinical database, T wave alternans levels are calculated in windows of 50 heartbeats as explained in section 3.6. In order to present results for the complete signals, the TWA levels in all windows are combined by calculating the 90th percentile of all values of windows containing TWA, resulting in one TWA value per signal. The detected levels of T wave alternans in the University Hospitals database are shown in figure 8 for both tensor decomposition methods.

**Figure 8.** Estimated levels of T wave alternans using PARAFAC2 and CPD in TWA group vs control group. The CPD results contain several outliers which are represented in red. Their value was capped at 60 μV (indicated with a dotted line) to increase clarity of the plotted results. (a) Estimated level of TWA using PARAFAC2 ( $p = 0.008$ ). (b) Estimated TWA level using CPD ( $p = 0.3768$ ).
Download figure:
Standard image High-resolution image

**Figure 8.** Estimated levels of T wave alternans using PARAFAC2 and CPD in TWA group vs control group. The CPD results contain several outliers which are represented in red. Their value was capped at 60 μV (indicated with a dotted line) to increase clarity of the plotted results. (a) Estimated level of TWA using PARAFAC2 ( $p = 0.008$ ). (b) Estimated TWA level using CPD ( $p = 0.3768$ ).
Download figure:
Standard image High-resolution image

The results are split according to the patient group (TWA vs control). When comparing between groups, both PARAFAC2 and CPD show higher TWA amplitudes for the TWA group than for the control group, but the differences are only significant with PARAFAC2 ( $p = 0.008$ ). Additionally we also analysed the number of windows which contain TWA in both groups. The number is higher for the TWA group for both methods, but the result is again only significant for PARAFAC2 ( $p = 0.03$ ). On average TWA is detected in double the windows of a signal in the TWA group than in the control group. Since TWA is linked to increased heart rates we also performed an analysis of the heart rate. This did however not reveal significant differences between both groups (mean heart rate of 96.8 ± 9 bpm in the TWA group versus 96.1 ± 13 bpm in the control group). Other signal characteristics such as signal length and number of channels did not differ significantly between both groups either.

To provide a suggestion for a clinically useful TWA threshold, we have used the values obtained by both CPD and PARAFAC2 to classify the signals as either containing TWA or not. The ROC curve of these classifications is shown in figure 9. The AUC values of CPD and PARAFAC2 are respectively 0.64 and 0.88. Further analysis of the results for PARAFAC2 shows that for a threshold of 13 μV the method obtains a sensitivity of 80 $\%$ and a specificity of 83 $\%$ , indicating that this might be a good cut-off for further use.

**Figure 9.** ROC curve for classification using CPD ( ${\rm AUC} = 0.64$ ) and PARAFAC2 ( ${\rm AUC} = 0.88$ ).
Download figure:
Standard image High-resolution image

**Figure 9.** ROC curve for classification using CPD ( ${\rm AUC} = 0.64$ ) and PARAFAC2 ( ${\rm AUC} = 0.88$ ).
Download figure:
Standard image High-resolution image

5. Discussion and conclusion

Tensors are a novel concept in ECG processing. In this paper, we present a method to detect T wave alternans using tensor decomposition methods. When comparing the two proposed tensor decompositions, CPD and PARAFAC2, we can conclude that PARAFAC2 has clear advantages over CPD. While CPD obtains good results in ideal circumstances (e.g. no noise and perfectly aligned T waves), these conditions are rarely met in real-life situations. This is also demonstrated by the superior results of PARAFAC2 on the non-artificial signals in the Physionet database and the clinical dataset from the University Hospitals Leuven. CPD fails to deliver good results in these circumstances, as is shown by the non-significant results on the patient data from the University Hospitals Leuven. CPD seems however slightly more robust to noise, based on figures 6(c) and (d). Its inability to deal with changing heart rates and movement makes it nevertheless unacceptable for real-life applications. While it is possible to control changing heart rates up to a certain degree, the results in figure 6(b) show that CPD results deteriorate quickly and become unreliable even with a minor change. Since small heart rate differences are unlikely to be avoided in real-life circumstances (even with control for heart rate), CPD would not be a good choice for use in clinical practice.

PARAFAC2 on the other hand yields good results both on artificial and clinical ECG signals. The difference is especially large in the presence of a T wave shift such as in figure 5(h). Since this occurs frequently in clinical signals, where a changing heart rate is expected, it is important that methods can deal with this type of variation. Figure 6(d) shows that the results of PARAFAC2 are less robust under high noise conditions. This indicates that it is presumably not an optimal method for ambulatory measurements. However, since many TWA tests occur in clinical environments where noise can be avoided to a certain extent, this is not likely to be a large obstacle in real life. The results shown on figures 8(a) and 9 confirm this. Special attention for the choice of segmentation window is however needed when dealing with patients that have conditions that can alter the T wave interval and timing, such as Long QT Syndrome. The results on the Physionet database seem to indicate that in extreme cases the proposed segmentation window is not sufficient to detect TWA. It would perhaps be a possibility to increase the length of the segmentation window in patients with known LQTS.

When comparing the results of the Physionet database with other methods, the tensor-based approach performs better than the other methods found in literature and obtains comparable results with the best-scoring method of the CinC challenge (Sieed and Hasan 2008) (0.87 versus 0.911). The Kendall coefficient only compares the ranking of the different signals, and gives no information about whether the estimation of the TWA magnitude is correct. For this a fully annotated standardized TWA database with correct estimations of TWA magnitude is necessary which does not exist at this moment. When analyzing the results of the clinical dataset from the University Hospitals Leuven, the difference between TWA positive signals and controls is statistically significant for PARAFAC2 ( $p = 0.008$ ). Note however that also in the control group a certain amount of TWA is still detected. This could be explained by a suboptimal choice of detection threshold (e.g. the minimal number of consequent sign changes for TWA presence). Full optimization of this parameter (using a database where the exact amount of TWA is known as explained earlier) could improve results here. Another possible explanation however is that the clinical spectral based TWA test does not detect very low levels of T wave alternans ( $<$ 10 μV). The spectral based method commercialized by Cambridge Heart was chosen as reference for the clinical dataset because it is one of the 2 official FDA-cleared methodologies and has extensively been studied (Verrier et al 2011). To verify such hypotheses the artificial signals containing fixed amounts of TWA could be analyzed by the clinical methods to compare the results. This could be an interesting approach for a future study. The suggested clinical threshold of 13 μV could also be validated in this way. A limitation of this study is the small size of the clinical database of only nine patients. In order to fully validate the clinical threshold and to verify whether it has prognostic utility, a larger prospective dataset is needed. This would also enable us to examine if the detected estimations of TWA magnitude correlate with clinical findings such as occurrences of arrhythmic events. A dedicated study would furthermore allow us to study the link between TWA levels and heart rate, as multiple studies have shown that TWA increases with increasing heart rate (Kavesh et al 1998, Kitamura et al 2002, Burattini et al 2014). The limited sample size did not allow us to make strong observations in this sense here.

The obtained results show that tensor methods are a promising tool in ECG signal processing. They are fully automated and require only a minimal amount of preprocessing. One disadvantage is the need for multilead ECG signals to construct the initial tensor since ECG is often measured using a single channel. In those cases the proposed methods can not be applied for TWA detection unless a different tensor construction approach is used. The presented method is fairly general and can (with slight adaptations) also be utilized in different applications. Other types of alternans such as ST alternans or QRS alternans can be evaluated easily by constructing a tensor that consists of other ECG segments. Another example is the detection of irregular heartbeats, of which the first results were presented already (Goovaerts et al 2015a).

Based on the obtained results we can conclude that the proposed methods are a robust way of detecting T wave alternans in different conditions. PARAFAC2 is a better choice as a tensor decomposition method than CPD since the more general structure supports time shifts in subsequent T waves.

Acknowledgments

Research supported by Bijzonder Onderzoeksfonds KU Leuven (BOF): Center of Excellence (CoE) #: PFV/10/002 (OPTEC); SPARKLE—Sensor-based Platform for the Accurate and Remote monitoring of Kinematics Linked to E-health #: IDO-13-0358; The effect of perinatal stress on the later outcome in preterm babies #: C24/15/036; TARGID - Development of a novel diagnostic medical device to assess gastric motility #: C32-16-00364. Fonds voor Wetenschappelijk Onderzoek-Vlaanderen (FWO): Project # G.0A55.13N (Deep brain stimulation). Agentschap Innoveren & Ondernemen (VLAIO): STW 150466 OSA +; O& O HBC 2016 0184 eWatch. iMinds Medical Information Technologies: Dotatie-Strategisch basis onderzoek (SBO- 2016); ICON: HBC.2016.0167 SeizeIT. Belgian Federal Science Policy Office: IUAP #P 7/19/ (DYSCO, ‘Dynamical systems, control and optimization’, 2012-2017). Belgian Foreign Affairs-Development Cooperation: VLIR UOS programs (2013-2019). EU: European Union’s Seventh Framework Programme (FP7/2007-2013): EU MC ITN TRANSACT 2012, # 316679; The HIP Trial: # 260777. Erasmus +: INGDIVS 2016-1-SE01-KA203-022114. European Research Council. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC Advanced Grant: BIOTENSORS (# 339804). This paper reflects only the authors’ views and the Union is not liable for any use that may be made of the contained information.

G G is a PhD fellow of IWT Vlaanderen.

R W has received research funding from Biotronik, Boston Scientific and Medtronic and speakers and consultancy fees from Biotronik, Boston Scientific, Medtronic, St Jude Medical and Sorin. R W is supported as postdoctoral clinical researcher by the Fund for Scientific Research Flanders.

Automatic detection of T wave alternans using tensor decompositions in multilead ECG signals

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

Abstract

1. Introduction