Computer Software tool for heart rate variability (HRV), T-wave alternans (TWA) and heart rate turbulence (HRT) analysis from ECGs

Background

Quantitative evaluation of cardiac risk prediction parameters requires implementation of complex computer algorithms. Various commercial applications provide cardiologists with tools that deliver diagnostic parameters, including heart rate variability (HRV) [1–3], heart rate turbulence (HRT) [4], T-wave alternans (TWA) [5–7] and many more. The above-mentioned commercial tools, however, often lack the ability to automatically compare the results with other parameters or to verify the applied methods, thus presenting doctors with a final result without being able to examine or adjust the employed algorithms, and thereby precluding further studies. In order to perform diagnoses based on various parameters, physicians often must use different, incompatible medical equipment and software tools. Some of these tools require advanced programming skills such as C-language programming to take advantage of the offered open source procedures. As a staring point for description of the program reported in this work we provide a short overview of the available toolkits for HRV, HRT and TWA [4] analyses.

The heart rate variability (HRV) toolkits are listed in Table 1 along with their software implementations and main features. All listed packages offer both time and spectral analyses of HRV, but only some of them implement methods derived from non-linear dynamic techniques. Only the Finnish Kubios HRV analysis software is equipped with a graphical user interface (GUI). The toolkit comes with an extensive user guide and provides user support via e-mail.

Heart rate turbulence (HRT) denotes transient fluctuations in the heart rate that directly precede and follow a ventricular premature beat (VPB). The term HRT was coined in 1999 by Schmidt et al, who defined the method in a seminal paper published in Lancet [13]. Currently, many clinical studies have confirmed the significance of HRT as an important predictor of sudden cardiac death risk in patients with myocardial infarction. The computations involved in measurement of HRT parameters are quite simple and rely on measuring geometric relations of the HR signal in the vicinity of the VPB beat. Regrettably, the web page www.h-r-t.org with free software has been discontinued. The HRT analysis program was patented in 2002. Recently, guidelines for HRT measurement standards and physiologic interpretation in clinical studies were published under the supervision of the International Society for Holter and Noninvasive Electrophysiology (ISHNE) [14].

In the case of T-wave alternans (TWA), to the authors’ best knowledge, there is only one TWA analyser toolkit made available as open source software. The toolkit that can be accessed from the Physionet web-page [15] is a suite of Matlab procedures for T-wave analysis [16], which implement the spectral method (SM) and the modified moving average method (MMA). A possible explanation for such limited availability of TWA analysis software tools is the fact that this cardiologic technique requires special ECG measurement protocols and apparatus to detect microvolt T-wave alternans. However, a considerable number of TWA analysis methods have been proposed. For an excellent and comprehensive review of 12 methods see Martinez and Olmos [17]. In 2010 the U.S. Food and Drug Administration gave market clearance for the microvolt TWA spectral method as a viable measure for risk of sudden cardiac death (SCD).

The indicated indices that proved important in risk stratification of patients (e.g., post- infarction patients), unfortunately, are not computed in a single software environment or tool. This makes an overall evaluation of a patient’s cardiac state a time consuming task. It is also difficult to find the relationship between the parameters since the examination registrations are not simultaneous and the conditions may vary significantly.

In this paper a computer program featuring a unified graphical user interface (GUI) for analysis of HRV, HRT and TWA is presented. All of the parameters are computed from the ECG recordings or, alternatively in the case of HRV, from previously computed RR series. At its current development stage, the software is designed predominantly for the purpose of scientific research rather than clinical diagnosis. Different file formats of ECG recordings can be loaded into the program. The user is presented with a graphical interface that enables plotting of ECGs and visual verification of intermediate steps of the analyses.

Material and Methods

TWA:

T-wave alternans is a phenomenon of beat-to-beat alternation of T-wave amplitude and duration occurring at an increased heart rate (up to 110 bpm). Numerous studies have established that TWA analysis is a useful tool for prediction of arrhythmias and SCD risk stratification. The mechanisms underlying alternations in repolarisation are highly dependent on pathophysiologic conditions [5,6]. The dominant mechanism is dispersion of depolarisation, which can be induced either functionally during rapid heart pacing or when metabolic impairment of ion pumps causes regional changes in ionic gradients as a result of acute myocardial ischemia. The clinically significant TWA values range from 2 to 10 μV [6]. Therefore, being invisible to the unaided eye, microvolt TWA must be analyzed using computer algorithms.

In the designed program, TWA analysis is performed both in frequency [5] and in time domain by applying Poincare maps [7]. The T-wave period for analysis can be selected by the user. The selection is marked graphically on the plot of ECG cycles superposed and aligned to previously determined QRS complexes. The amplitude of the T-wave is measured for 7 equally spaced time instances – the first corresponding to the start, and the last to the end of the T-wave selected by the user as shown in Figure 2.

In the next step, for each of the seven T-wave phases, a series x_n of successive samples is constructed. In the frequency domain analysis method, each of the series undergoes the Fourier Transform. The level of the amplitude spectrum at frequency equal to 0.5 cycles per beat indicates the level of beat-to-beat alternation of the T-wave amplitude at a given phase. The alternans voltage and ratio indices are calculated according to the formulae 1 and 2, respectively:

where:

For the time domain analysis the scatter graph of the Poincare mapping is applied. The scatter graph is a plot of x′n+1 versus x′n, where x′ is the first order differential of series xn. In case of ECG registration of a patient unaffected by TWA, the points on the plot form an unstructured cluster. Otherwise, the points are situated in the second and fourth quarter of the scatter graph (Figure 3B). The DTWA index, defined in equation 3, can be calculated as the distance between the two clusters formed by even and odd scatter graph point pairs.

Values of both spectral V_alt and temporal D_TWA are computed for all six precordial leads.

By using the Poincaré mapping as well as spectral analysis methods, we have constructed spatio-temporal imaging of TWA content across V1–V6 precordial leads for m=7 T-wave sections. First, a matrix of size 7×6 DTWA indices is obtained. Figure 4 explains the procedure for construction of this matrix. Then, in order to generate a smooth image of such a spatio-temporal distribution of TWA, a biharmonic spline interpolation technique is employed for spanning an interpolation surface supported by 42 knots at space-time locations of the precomputed DTWA indices. This type of imaging offers a new qualitative way of viewing TWAs. Example spatio-temporal maps for TWA+ patients are shown in Figure 3C,D. Hot colors indicate regions of high TWA levels that can be used for both spatial (in terms of precordial lead location) and temporal (T-wave section number) identification of heart repolarization abnormalities. The program provides the possibility to analyse both sinus rhythm and paced rhythm recordings.

HRT:

Heart rate turbulence is a response of the sinus node to the ectopic beat, characterized by rapid acceleration prior to ventricular premature beat (VPB), followed by deceleration of the heart rate. The most probable mechanism underlying HRT is the baro-reflex autoregulation process, which responds to the deviation in blood pressure caused by ectopic beat occurrence [13]. Missing or weakened response of such a nature reflects a malfunctioning of the autonomic control system and is considered a powerful risk indicator in patients after myocardial infarction.

In the designed program, HRT onset and HRT slope parameters are computed using the Schmidt algorithm [13]. Turbulence onset (TO) is defined (eq. 4) as the mean percentage difference between the average interbeat interval (IBI) after and before a ventricular premature beat.

where E[.] denotes the expectation operator and n is the index of the pre-extrasystole interval.

The turbulence slope (TS) is the maximum slope of the linear regression line for each sequence of 5 consecutive normal intervals taken from the 15 IBIs following the VPB, providing they are free from arrhythmia and artifacts.

Because the HRT pattern can be masked by heart rate variability, these parameters are usually calculated as an averaging responses to >5 ventricular premature contractions. Long-term Holter recording is recommended to assess this parameter [14].

HRV:

Heart rate variability (HRV) is a measure of beat-to-beat variations in the length of the RR interval, and is a marker of autonomic nervous system modulation of the sinus node activity. Limited studies have linked short-term HRV to sudden death, and its use for risk stratification is currently not recommended [24].

The designed program enables computation of statistical parameters in time domain (mNN, SDNN, RMSSD, pNN50) and a number of parameters in the frequency domain (LF and HF power, LF/HR ratio). The time domain parameters are calculated directly from the RR series. The mNN and SDNN are mean and standard deviation of normal-to-normal RR intervals, RMSSD is the root mean square of the differences between consecutive RR intervals, and pNN50 is the percentage of a number of consecutive RR intervals differing by more than 50 ms. Spectral analysis based on the Fourier Transform (FT) requires preprocessing of the RR interval sequence. Artefacts, including ectopic beats [25], and missed and false detections, are located in the series by applying thresholding to the second derivative of the HRV series and then removed. In order to obtain uniformly sampled series, the RR samples are interpolated and the result is re-sampled with an appropriate frequency. A frequency of 2Hz is assumed as the maximum possible frequency of occurrence of the heart beat event. Interpolation is performed using the cubic spline method, which is based on piecewise polynomial functions of the third order. The frequency characteristic is determined using the parametric autoregressive moving average (ARMA) parametric model and nonparametric (periodogram) model, (i.e. the squared value of the Fourier Transform of the analysed HRV) [26]. In the implemented Welch periodogram method, the signal is divided into a number of overlapping intervals and multiplied by a window function prior to FT analysis. The periodograms are computed for each interval and averaged. By applying such a procedure, the obtained spectrum is smoother than in the case of classic periodogram. By adjusting the degree of overlap, the number of intervals, and type of the window, it is possible to improve the resolution of the spectral characteristic at the cost of its smoothness [26].

In non-parametric methods the frequency characteristics are computed directly from the examined signal. In parametric methods, on the other hand, we are looking for a hypothetical system that generates the examined signal. Since the heart’s natural pacemaker follows complex periodic dynamics, many simplifications need to be introduced into the model. Nevertheless, a good approximation of the underlying mechanisms controlling heart rate variability can be achieved with the use of such models.

In the moving average modeling (MA) approach, the HRV series is modelled by a weighted linear combination of s(k), k=1,…,q samples of a Gaussian noise process. The modeling task is to determine the vector of weight coefficients b_k that would generate the best fit to the analyzed HRV series y(n) according to the formula:

One of the model design parameters is the order of the model given by q. Too low an order of the model would omit important HRV details, whereas too high an order can yield a noisy series that would be difficult for meaningful interpretation.

In the autoregressive modeling (AR) approach, the HRV series is modeled by a recursive process in which the current model sample is generated from a predetermined number of its previous samples (hence the model name “autoregressive”) multiplied by appropriate coefficients ak, k=1, …, p. The corresponding equation is the following:

Again, the model order p in an important parameter that has to be decided in order to obtain a suitable AR model of the HRV series.

Both described models have limitations. The MA model is driven by Gaussian noise and omits inherent dynamics of the modelled object. On the other hand, the AR model does not take into account any stochastic mechanisms that influence the HRV. Hence, a modeling approach that combines the two models is a frequently adopted technique in HRV parametric modeling. The resulting model is termed the autoregressive moving average model (ARMA) and the resulting model equation assumes the form:

Here, the goal of the modelling process is to find vectors of coefficients ak and bk so that the output y(n) of the system would be as close as possible to the examined HRV series. This is a challenging task because in practice 10–20 parameters need to be simultaneously determined. There have been, however, special algorithms developed to solve this multicriterion optimisation problem. The most effective ones are the Burg’s and Durbin’s algorithms [27].

Another, indicated problem is the selection of the orders p,q of the model. These orders influence the complexity of the system generating the HRV series. Assuming the same orders for all the patients or even for the same patient at various conditions should be considered too far a simplification. The algorithms implemented in the presented software provide means for automatic computation of the model order by optimising model complexity versus the modelling error. The implemented order selection algorithm is based on the so-called Akaike’s information criterion, which indicates the goodness of fit of the model to the experimentally measured HRV series [27].

The interested reader may probe further into AR, MA and ARMA modelling by studying classic texts on the subject [26, 28].

Results

A user-friendly software package for quantitative evaluation of HRV, HRT and TWA in ECGs was developed. The software tool offers the following main features for ECG processing and analysis:

ECG preprocessing tools and file handling tools:

ECG analysis tools:

Clinical recordings that have already been examined by this software package were provided by the Department of Electrocardiology of the Medical University in Lodz. Programs used in this software package were used to assess TWA in postinfarction patients with low ejection fraction (EF%) and preserved left ventricular function [29], HRV in postinfarction patients with impaired left ventricular function [30], and HRV in patients with chronic heart failure and implanted cardiac resynchronization devices [31]. Another group of recordings come from the Physionet online ECG signal repository [15].

Each of the supported analyses can be performed using different methods with various parameters, which makes it possible to cross-validate the results obtained. The derived parameters (given in time, frequency and spatio-temporal domain) are visualized for easier interpretation. Preview of the intermediate steps of the analysis (QRS detection, HRV artefact correction, etc.) is a good and easily available source of information on possible errors such as beat misdetections. Thanks to this previewing functionality, the parts of registrations prone to errors (eg, fragments with very frequent occurrence of ectopic beats) can be excluded from further analyses.

Example results of a TWA analysis are shown in Figure 3. Spectra of the series of T-wave amplitudes for recordings from the precordial leads are plotted in Figure 3A. Results of the temporal analysis obtained by applying the Poincare maps are shown in Figure 3B–D illustrate the interpolation maps of TWA voltage obtained from spectral and spatio-temporal analysis, respectively. Similar results can be obtained using the two methods, while the computational cost of the latter is much lower [7]. To achieve interpretable results of the spectral analysis, a relatively long analysis period is required.

In spectral analysis of HRV the choice of the method greatly influences the results. Autoregressive modelling is suitable for estimating spectra containing sharp peaks. Modelling the bands of the smooth spectra is inefficient with the AR method. On the other hand, the moving average (MA) method is more precise in modelling relatively smooth spectra. ARMA modelling being the combination of both autoregressive and moving average models is capable of estimating spectra partially rich with peaks and partially smooth. The cardiologist analyzing the results of the HRV modelling should be aware of the properties of these algorithms. The shape of the HRV spectra is only partially dependent on the nature of autonomous nervous system behaviour; it also depends on the applied model type and model order (Figure 9). In the designed program, an automatic model order selection algorithm is implemented and the results of different methods are available for comparison and cross-validation.

Discussion

Although isolated studies have demonstrated the utility of many different variables for sudden cardiac death (SCD) risk stratification (imaging-based, ECG-based, exercise test-based, Holter-based and others), the ideal tool to identify patients who might experience life- threatening arrhythmias has not been found. Among the risk factors assessed non-invasively, low left ventricular ejection fraction (LVEF%) is the most powerful predictor of serious events in ischemic and non-ischemic cardiomyopathy. In fact, it is the sole variable determining the recommendation for implantable cardioverter-defibrillator (ICD). However, the majority of SCDs occur in patients with more preserved or normal EF, and the evaluation of risk in this large group of patients is justified by means of the secondary factors computed within the proposed software tool [24]. Assessment of risk factors derived from the most frequently performed examinations, such as supine and exercise ECG, can help to resolve the problem of risk stratification for SCD in patients without left ventricular impairment.

Conclusions

The developed software package is a flexible and comprehensive tool for TWA, HRT and HRV quantitative analysis from ECGs. It offers the possibility to choose from various modelling methods and to modify parameters of the applied algorithms. The interface is user-friendly and enables interactive selection of ECG signal fragments for analysis. Each software function has been developed and tested in close cooperation with cardiologists. The authors believe that the proposed software tool for ECG analyses can offer a basis for a better common understanding of engineers and clinicians working together in the difficult field of cardiology.