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

BACKGROUND: This paper presents a software package for quantitative evaluation of heart rate variability (HRV), heart rate turbulence (HRT), and T-wave alternans (TWA) from ECG recordings. The software has been developed for the purpose of scientific research rather than clinical diagnosis. MATERIAL/
METHODS: The software is written in Matlab Mathematical Language. Procedures for evaluation of HRV, HRT and TWA were implemented. HRV analysis was carried out by applying statistical and spectral parametric and nonparametric methods. HRT parameters were derived using the Schmidt algorithm. TWA analysis was performed both in spectral and in time domain by applying Poincare mapping. A flexibility of choosing from a number of classical modelling approaches and their modifications was foreseen and implemented. The software underwent preliminary verification tests both on ECGs from the Physionet online ECG signal repository and recordings taken at the Department of Electrocardiology of the Medical University Hospital in Lodz.
RESULTS: The result of the research is a program enabling simultaneous analysis of a number of parameters computed from ECG recordings with the use of the indicated analysis methods. The program offers options to preview the intermediate results and to alter the preprocessing steps.
CONCLUSIONS: By offering the possibility to cross-validate the results of analyses obtained by several methods and to preview the intermediate analysis steps, the program can serve as a helpful aid for clinicians in comprehensive research studies. The software tool can also be utilized in training programs for students and medical personnel.

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.

Table 1

Software tools for HRV analyses****.

Ref	Toolkit	Platforms	Time analysis	Spectral analysis	Preprocessing or other analyses	GUI
[8]	J.E. Metus, A.L.Golderberg, “Heart Rate Variability Analysis With The HRV Toolkit” www.physionet.org/tutorials/hrv-toolkit	C-language, (UNIX, Windows, Mac Os, Solaris)	AVNN, SDNN, SDANN, SDNNIDX, rMSSD, pNN50	TOTPWR, ULF, VLF, LF, HF, LF/HF	Tool for detection and elimination of outliers	−
[9]	J. Niskanen et. al, “Kubios-HRV” http://kubios.uku.fi/	Matlab standalone application (Windows, Linux)	SDNN, SDSD, RMSSD, pNN50, TINN	AR estimates of VLF, LF, HF, LF/HF	Poincare plots, ApEn, SampEn, DFA, CorrDim	+
[10]	M. Lado, A. Mendez, L. Rodriguez-Linares, X. Vila, “R-HRV” www.cran.r-project.org/web/packages/RHRV/	C-language (Mac OS X, Windows)	−	VLF, LF, HF, LF/HF	Outlier filtering, interpolation	−
[11]	D. Kaplan, P. Staffin “Software For Heart Rate Variability” www.macalester.edu/~kaplan/hrv/doc/	Matlab m-files	SDNN, RMSSD, pNN50, TINN	FFT estimates of VLF, LF, HF, LF/HF	Outlier handling, ApEn*, DFA**	−
[12]	Nevrokard, “aHRV”	Windows (free demo)	AVNN, SDNN, SDANN, SDSD, RMSSD, NN50,	AR and FFT spectrum	Poincare plots, CZF*** Analysis	+

ApEn – Approximate entropy;

DFA – Detrended Fluctuation Analysis;

CZF – Conte, Zbilut, Federici Analysis;

see section on HRV for explanation of other symbols.

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 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

The software is written in Matlab Mathematical Language. The test version of the program can be downloaded without charge from . In order to run the program it is necessary to install Matlab Compiler Runtime (MCR) first. A suitable guide is provided for installing the program. The graphical user interface of the program is shown in Figure 1.

Figure 1

Graphical User Interface of the designed program.

The recorded ECG signals or pre-computed RR series can be loaded from binary, text or system-specific format files (e.g. digital ECG Medea, Gliwice, Poland [18], and Finometer, Finapres Medical Systems, Amsterdam, the Netherlands [19]). Suitable filtering procedures are available for denoising the recorded ECGs and beat-by-beat synchronous averaging.

Baseline wandering removal is implemented by detecting signal samples corresponding to PQ segments followed by cubic spline interpolation [20-22]. The QRS complexes are detected using digital filters based on the Haar wavelet [23]. A Haar filter of 20th order is a bandpass linear-phase filter whose passband is centered at f=18.5 Hz (for sampling frequency fs=500 Hz), while its first and second zeros correspond to 0 Hz and 50 Hz, respectively. Therefore, such a filter retains signal frequency components of QRS complexes and suppresses other signal components such as base-line wander and 50 Hz interference coming from the electric powerline network. The reference PQ-interval points are found 66msec prior to the Q-wave’s maximum downslope. The software enables plotting ECG records, visualizing QRS detections, computing and plotting HRV series.

The current version of the software enables evaluation of parameters that are claimed to be important for SCD risk stratification, mortality after myocardial infarction and for assessment of autonomic system dysfunction in patients with diabetic neuropathy or hypertension. However, further studies are required to evaluate the true validity of these parameters in clinical practice [24].

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.

Figure 2

T-wave range selection.

In the next step, for each of the seven T-wave phases, a series x 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:

V –alternans voltage,

k – alternans ratio,

S0.5 – level of spectrum at frequency 0.5 cycles per beat,

Snb and σnb– mean level and standard deviation of spectrum at the noise band defined as 0.44–0.49 cycles per beat.

For the time domain analysis the scatter graph of the Poincare mapping is applied. The scatter graph is a plot of x′+1 versus x′, where x′ is the first order differential of series x. 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 D index, defined in equation 3, can be calculated as the distance between the two clusters formed by even and odd scatter graph point pairs.

Figure 3

Graphical representation of the results of TWA analysis: Spectrum (A) and Poincare map (B) for an arbitrarily chosen lead and section; interpolation maps based on spectral (C) and spatiotemporal (D) parameters.

Values of both spectral V and temporal D 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 D 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 D 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.

Figure 4

Construction of a spatio-temporal image of TWA for precordial ECG leads by means of interpolation (the interpolation knots are indicated by gray circles).

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 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 a and b 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:

Loading ECG recordings from text or binary files.

Plotting the loaded ECG recordings of every lead as well as the result of QRS detections accompanied by the preview of RR series. The preview plots of the currently loaded ECGs are built-in within the program’s GUI (Figure 1). The plots can also be opened in a new window.

Detection of QRS complexes of the SAN (sinoatrial node) rhythm as well as the rhythm stimulated by artificial pacemakers (e.g. atrial sensed, ventricular paced rhythm) (Figure 5).

Figure 5

Paced rhythms QRS detection

Choosing either the entire ECG recording or its fragment for further analyses.

Denoising ECGs with synchronized averaging method as shown in Figure 8 (averaging a number of ECG cycles synchronized by the QRS complex). Such an approach enables reduction of random noise signal components and filtering out ECG characteristic features repeating at every averaged cycle (e.g. hidden P wave).

Figure 8

Synchronized averaging.

Exporting the loaded signals (either source ECG recordings or preprocessed ECGs) to .txt and .mat files.

ECG analysis tools:

HRT analysis. The recommended number of 15 artefact-free cycles following VPB can be selected by the user. The values of TO and TS parameters are accompanied by the plots presenting the algorithm’s detection of VPB beats fulfilling the conditions described in the section about HRT analysis. The user is also provided with the ECG plot where the mentioned detections are marked. This allows a more comprehensive and reliable analysis. A typical outcome of HRT analysis is shown in Figure 6.

Figure 6

HRT analysis.

TWA analysis (Figure 3). After graphical selection of T-wave start and end points (Figure 2) the user is provided with both spectral and Poincare analysis plots for each of the precordial leads at every time section. Finally, the spatio-temporal maps described in previous sections are displayed.

HRV analysis. As shown in Figure 7, the user can change the preprocessing options. The sensitivity of artefact correction can be set by the user, which may be found useful in case of highly contaminated recordings or in situations were certain natural rhythm changes are wrongly interpreted by the algorithm as artefacts. Statistical HRV parameters are output to the console. The parameters, including type of analysis, model type, and model order are also user-defined.

Figure 7

HRV analysis.

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.

Figure 9

Spectra of the same HRV series obtained by parametric modelling with different model types and orders.

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.