Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3523127

Formats

Article sections

- Abstract
- I. Introduction
- II. Algorithmic Framework
- III. Validation and Performance Assessment
- IV. Discussion and Conclusions
- References

Authors

Related links

IEEE Trans Biomed Eng. Author manuscript; available in PMC 2013 October 1.

Published in final edited form as:

Published online 2012 August 2. doi: 10.1109/TBME.2012.2211356

PMCID: PMC3523127

NIHMSID: NIHMS419572

Dept. Anesthesia, Massachusetts General Hospital – Harvard Medical School, Boston, MA, USA, and with the Dept. Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA

Luca Citi: ude.tim.tatsoruen@iticl; Emery N Brown: ude.tim.tatsoruen@bne; Riccardo Barbieri: ude.tim.tatsoruen@ireibrab

The publisher's final edited version of this article is available at IEEE Trans Biomed Eng

See other articles in PMC that cite the published article.

The presence of recurring arrhythmic events (also known as cardiac dysrhythmia or irregular heartbeats), as well as erroneous beat detection due to low signal quality, significantly affect estimation of both time and frequency domain indices of heart rate variability (HRV). A reliable, real-time classification and correction of ECG-derived heartbeats is a necessary prerequisite for an accurate on-line monitoring of HRV and cardiovascular control. We have developed a novel point process based method for real-time R-R interval error detection and correction. Given an R-wave event, we assume that the length of the next R-R interval follows a physiologically motivated, time-varying inverse Gaussian probability distribution. We then devise an instantaneous automated detection and correction procedure for erroneous and arrhythmic beats by using the information on the probability of occurrence of the observed beat provided by the model. We test our algorithm over two datasets from the Physionet archive. The Fantasia normal rhythm database is artificially corrupted with known erroneous beats to test both the detection and correction procedure. The benchmark MIT-BIH Arrhythmia database is further considered to test the detection procedure of real arrhythmic events and compare it with results from previously published algorithms. Our automated algorithm represents an improvement over previous procedures, with best specificity for detection of correct beats, as well as highest sensitivity to missed and extra beats, artificially misplaced beats, and for real arrhythmic events. A near-optimal heartbeat classification and correction, together with the ability to adapt to time-varying changes of heartbeat dynamics in an on-line fashion, may provide a solid base for building a more reliable real-time HRV monitoring device.

Heart rate variability (HRV) techniques [1]–[3] provide a window into the many physiological factors that modulate the normal heart rhythm. In particular, they have been found useful for non-invasive autonomic tone assessment in a wide range of clinical and non-clinical scenarios. In order to provide reliable results, these techniques require uninterrupted series of normal R-R intervals. Peak detection errors – when the algorithm misses a beat and/or detects one when there is none – and ectopic beats often determine abrupt changes in the R-R interval series that can lead to substantial deviations of the HRV indices.

The increasing quality of electrocardiogram (ECG) recording systems together with better R-peak detection algorithms has reduced the problem of noisy R-R series in controlled clinical and research environment. Nevertheless, this issue is still of primary importance in new applications of HRV studies, such as ambulatory monitoring in very unstructured environments and extreme conditions, and when the ECG is acquired using single-channel miniature devices, smartphones, or wearable technology.

The simple approach of discarding all sections contaminated by irregular R-R intervals might pose some problems. For example, if the rate of occurrence of mis-detections or ectopic events is not uniform but it is dependent on specific parts of the protocol (e.g., those with frequent motion artifacts) or on the physiological state of the subject (e.g., states that may result in abnormal morphology of the QRS or in higher incidence of ectopic events), then the exclusion of these periods might skew the resulting HRV indices [4]–[6]. When the goal of the research is to study HRV in extreme conditions (e.g., during rescue operations, military missions, or hiking expeditions) removing all sections contaminated by R-peak misdetections caused by motion artifacts is simply not an option.

A better approach is to attempt to correct the corrupted R-R series in order to obtain a new artifact-free R-R series reflecting the underlying HRV dynamics. This new R-R series can then be used for further processing using the HRV analysis techniques of choice. To date, the detection and correction of irregular beats has been mainly achieved by direct human expert evaluation and, more recently, by automatic techniques [4]–[16].

We have developed a novel point process based method for real-time R-R interval error detection and correction. Given an R-wave event, we assume that the length of the next R-R interval follows a time-varying history-dependent inverse Gaussian (IG) probability distribution [17]. We then use this model for the detection and correction of erroneous and arrhythmic beats. A preliminary version of the detection step of this novel approach has been briefly sketched in our previous publication [18]. In Section II-C we present the detection algorithm in detail and later in the paper we show its performance in comparison with three established methods for the detection of erroneous (Section III-A) and ectopic (Section III-C) beats. Whenever a beat is classified as erroneous, the correction step attempts to fix it according to the point-process model. In Section II-D we present the correction algorithm and in Section III-B we assess its performance in comparison with four reference methods.

A software tool that implements the algorithm described in this manuscript can be tested from our website [19]. It can be used as a preprocessing stage to detect and correct erroneous and ectopic beats in R-R series before analysing it with the chosen HRV technique [1]–[3], such as time-frequency analysis (LF/HF, …), standard time-domain methods (SDNN, SDANN, RMSSD, …), and non-linear indices (DFA, SampEn, MSE, …). Some methods, such as the Lomb spectrogram [20] do not require the correction step but might still benefit from the accurate detection capabilities of our algorithm.

Each cardiac contraction is initiated by the synchronous depolarization of the heart’s pacemaker cells beginning in the sinoatrial (SA) node of the right atrium and then propagating through its specialized conduction system to the left atrium and to the two ventricles. Following every depolarization, the transmembrane potentials of these cells return to their resting potentials and begin anew their spontaneous rise toward threshold [21]. Deterministic models of this integrate (rise of transmembrane potential)-and-fire (depolarization) mechanism, such as the integral pulse frequency modulation (IPFM) model, are used regularly to simulate heart beats [22]–[24] and for the analysis of HRV [25].

An elementary, stochastic integrate-and-fire model is the Gaussian random walk model with drift. We assume that the afferences to the SA node can be represented by random walk excitatory inputs responsible for the basal cardiac rhythm (B), and random walk excitatory (E) and inhibitory (I) inputs, reflecting the influence of the autonomic nervous system through the sympathetic and para-sympathetic branches. In this case, the membrane voltage equation is governed by:

$$dV(t)={\alpha}_{\text{B}}\phantom{\rule{0.16667em}{0ex}}{dN}_{\text{B}}(t)+{\alpha}_{\text{E}}\phantom{\rule{0.16667em}{0ex}}{dN}_{\text{E}}(t)-{\alpha}_{\text{I}}\phantom{\rule{0.16667em}{0ex}}{dN}_{\text{I}}(t)$$

where *N*_{B}(*t*) ~ *P* (*λ*_{B}), *N*_{E}(*t*) ~ *P* (*λ*_{E}) and *N*_{I}(*t*) ~ *P* (*λ*_{I}) are independent Poisson processes that govern the respective times of the basal, excitatory and inhibitory inputs, while *α*_{B}, *α*_{E} and *α*_{I} are the magnitude of steps up or down [26], [27]. Assuming *V* (0) = 0, *m*_{0} = *α*_{B}*λ*_{B}, *m*_{1} = *α*_{E}*λ*_{E} − *α*_{I}*λ*_{I}, and
${\sigma}^{2}={\alpha}_{\text{B}}^{2}{\lambda}_{\text{B}}+{\alpha}_{\text{E}}^{2}{\lambda}_{\text{E}}+{\alpha}_{\text{I}}^{2}{\lambda}_{\text{I}}$, and using the properties of Poisson processes, the mean of *V* (*t*) is E [*V* (*t*)] = (*m*_{0} + *m*_{1}) *t* [22], [23] and its variance Var[*V* (*t*)] = *σ*^{2}*t*. A diffusion approximation of *dV* (*t*) before hitting the threshold *θ* is given by the Wiener process with drift:

$$dV(t)=({m}_{0}+{m}_{1})\phantom{\rule{0.16667em}{0ex}}dt+\sigma \phantom{\rule{0.16667em}{0ex}}dW(t)$$

where *W*(*t*) is a standard Wiener process (Brownian Motion). The probability density of the first passage time for this process, i.e., the times between threshold crossings, is given by the inverse Gaussian probability density [28], [29] defined as:

$$f(t)=\frac{\theta}{\sqrt{2\pi {\sigma}^{2}{t}^{3}}}{\text{e}}^{-\frac{{(\theta -({m}_{0}+{m}_{1})t)}^{2}}{2{\sigma}^{2}t}}$$

which can be written in the usual form after performing the change of variables *μ* = *θ*/(*m*_{0} + *m*_{1}) and *λ* = *θ*^{2}/*σ*^{2}. The inverse Gaussian probability density has been used to study both neural processes [30] and heart beats [17], [31].

An ordered set
${\{{u}_{j}\}}_{j=1}^{J}$ of consecutive heartbeat times (for example, but not limited to, R-waves detected from an ECG) recorded in an observation interval (0, *T*] can be interpreted as a realization of a point process [32], a type of stochastic process for which any one realization consists of a set of isolated points. For *t* (0, *T*], the cádlág function *N*(*t*) = max{*k: u _{k}* ≤

In this paper we assume a history dependent, time-varying, model of the the probability distribution of the waiting time until the next R-wave according to a physiologically motivated (see Section II-A) inverse Gaussian distribution. Specifically, given any beat event *u _{k}*, the probability density function (PDF) of the length of the next R-R interval,

$$f(\tau -{u}_{k}\mid \mu \phantom{\rule{0.16667em}{0ex}}({\mathit{H}}_{k},\mathit{\theta}(t)),\lambda (\mathit{\theta}(t)))=\sqrt{\frac{\lambda (\mathit{\theta}(t))}{2\pi {(\tau -{u}_{k})}^{3}}}exp\left\{-\frac{1}{2}\frac{\lambda (\mathit{\theta}(t))\phantom{\rule{0.16667em}{0ex}}{(\tau -{u}_{k}-\mu \phantom{\rule{0.16667em}{0ex}}({\mathit{H}}_{k},\mathit{\theta}(t)))}^{2}}{(\tau -{u}_{k})\phantom{\rule{0.16667em}{0ex}}{\mu}^{2}}\right\}$$

where the shape parameter, *λ*, and the mean of the distribution, *μ*, depend on a vector of time-varying parameters ** θ**(

$$\mu \phantom{\rule{0.16667em}{0ex}}({\mathit{H}}_{k},\mathit{\theta}(t))=\sum _{i=1}^{P}{\theta}_{i}(t)\phantom{\rule{0.16667em}{0ex}}{w}_{k-i}.$$

(1)

A local maximum likelihood method [17] is used to estimate the unknown time-varying parameter set ** θ**(

$$L(\mathit{\theta}(t)\mid {U}_{m:n})=\sum _{k=m+P}^{n-1}\omega (t-{u}_{k+1})log\phantom{\rule{0.16667em}{0ex}}[f({u}_{k+1}-{u}_{k}\mid \mu ({\mathit{H}}_{k},\mathit{\theta}(t)),\lambda (\mathit{\theta}(t)))]$$

(2)

where *ω* (*τ*) = e^{−}* ^{ατ}* is an exponential weighting function for the local likelihood [17], [33].

After observing the (*k*+1)-th R-wave, at time *u _{k}*

A straightforward way to test the hypothesis that the beat at time *u _{k}*

One hypothesis that we consider is that the beat at time *u _{k}*

Another hypothesis is that between the events *u _{k}* and

$$\underset{{u}_{k}}{\overset{\tau}{\int}}f(\tau -{\tau}^{\prime}\mid {\mu}_{2}({\tau}^{\prime}-{u}_{k}),{\lambda}_{2})f({\tau}^{\prime}-{u}_{k}\mid {\mu}_{1},{\lambda}_{1})\phantom{\rule{0.16667em}{0ex}}\text{d}{\tau}^{\prime}.$$

(3)

Rather than attempting to solve (3), here we make the simplifying assumption that the resulting probability density function can be approximated by an IG distribution with mean *μ*_{1+2} = *μ*_{1} + *μ*_{2}(*μ*_{1}) and shape parameter

$${\lambda}_{1+2}={\lambda}_{1}\frac{{({\mu}_{1}+{\mu}_{2}({\mu}_{1}))}^{3}}{{(1+{\stackrel{\sim}{\theta}}_{1}({u}_{k}))}^{2}{\mu}_{1}^{3}+{\mu}_{2}^{3}({\mu}_{1})}$$

(4)

whose mathematical derivation is shown in the Appendix. Using these results, we can now evaluate the log probability density of observing a beat at *u _{k}*

A further hypothesis is that the beat at time *u _{k}*

This hypothesis accounts for the case that two beats in a row, *u _{k}*

The last hypothesis accounts for resetting ectopic beats, i.e., ectopic beats not followed by a compensatory pause. As a result of the absence of the compensatory pause, not only the R-R interval *u _{k}*

As the algorithm uses a sliding window of duration *W* to fit the model (see end of Section II-B), no model is available before time *W* from the beginning of the recording. A simpler bootstrap detection algorithm is used during this period to provide some results even for beats *u _{k}* (0,

After an outlying beat is detected, our algorithm attempts to correct it. The correction strategy depends on the type of erroneous beat detected, as reported in detail below. Finally, the attempted solution is accepted only if it is an improvement according to the improvement check that we describe in the following section.

The case of an extra beat is simply addressed by removing the beat.

In case a missed beat is detected, a new beat is inserted between the event *u _{k}* and the event

$${\tau}^{\prime}=\underset{{\tau}^{\prime}\in ({u}_{k},{u}_{k+1})}{argmax}f({\tau}^{\prime}-{u}_{k}\mid {\mu}_{1},{\lambda}_{1})\phantom{\rule{0.16667em}{0ex}}f({u}_{k+1}-{\tau}^{\prime}\mid {\mu}_{2}({\tau}^{\prime}-{u}_{k}),{\lambda}_{2})$$

(5)

where *μ*_{1}, *λ*_{1}, *μ*_{2}(*τ*′ − *u _{k}*) and

The case of a misplaced beat is similar to the previous case. The event at time *u _{k}*

$${\tau}^{\prime}=\underset{{\tau}^{\prime}\in ({u}_{k},{u}_{k+2})}{argmax}f({\tau}^{\prime}-{u}_{k}\mid {\mu}_{1},{\lambda}_{1})\phantom{\rule{0.16667em}{0ex}}f({u}_{k+2}-{\tau}^{\prime}\mid {\mu}_{2}({\tau}^{\prime}-{u}_{k}),{\lambda}_{2}).$$

(6)

Equations (5) and (6) are solved numerically using a Newton-Raphson algorythm.

To correct the case of two misplaced beats and find the optimal times of the two unknown beats we use a *greedy algorithm*. Each one of the beat times is optimized in turn, using (6) for the first beat and a similar equation for the second beat, while keeping the other fixed. This process is iterated until convergence.

The case of a resetting beat cannot be corrected like the other cases because it would require shifting all the following beats. Therefore our algorithm defers the decision about the action to take to the user, who can choose the best strategy for the type of HRV analysis that follows (for example omit the beat or suspend adaptive causal estimation).

In the cases above, the correction step changes the series of events *U* = {*u _{j}*} into an alternative series

Before the new series is accepted, we evaluate the probability of observing *Q* events following *u _{k}* for both series:

$$\begin{array}{l}\mathbf{P}({\widehat{U}}_{k-P;k+Q})=\sum _{j=k}^{k+Q-1}log\phantom{\rule{0.16667em}{0ex}}\left[f\left({\widehat{u}}_{j+1}-{\widehat{u}}_{j}\mid {\widehat{\mathit{H}}}_{j},\stackrel{\sim}{\mathit{\theta}}({u}_{k})\right)\right],\\ \mathbf{P}({U}_{k-P:k+Q})=\sum _{j=k}^{k+Q-1}log\phantom{\rule{0.16667em}{0ex}}\left[f\left({u}_{j+1}-{u}_{j}\mid {\mathit{H}}_{j},\stackrel{\sim}{\mathit{\theta}}({u}_{k})\right)\right].\end{array}$$

The corrected series is accepted only if **P**(*Û _{k}*

Empirical sensitivity analysis conducted on in-house recordings showed that the model provides robust results for parameters within meaningful ranges of values.

Without conducting an exhaustive research on the whole parameter space, we found that the following settings provide excellent results in most conditions: *P* = 5, *W* = 60 s, *α* = 0.02 s^{−1}, *η*_{e} = 3, _{e} = 8, *η*_{s} = 0, _{s} = 4, *η*_{m} = 2, _{m} = 7, *η*_{t} = 8, _{t} = 28, *η*_{r} = 6, _{r} = 14, and *Q* = 3. These values of the parameters were also used during the analysis described in the following of this manuscript.

We tested the performance of our detection algorithm on a set of real R-R timeseries with artificially deleted, shifted, and inserted beats and compared our results with three well-established approaches.

The first method that we implemented for comparison is the well-established algorithm by Berntson et al. [9] which uses two criteria to assess whether a heartbeat is an artifact: Maximum Expected Difference (MED) and Minimal Artifact Difference (MAD). In order to allow the algorithm to adapt to non-stationary data, for each R-event *u _{k}*, MED and MAD were estimated using events in a 60 s sliding window preceding

For this comparison, we used a subset of the Fantasia database [34] from the Physionet archive [35]. We selected all the records that, according to the annotation database, contained no more than two non-“N” heartbeats.^{1} We tested all four algorithms on each pristine record and recorded the number of “N” heartbeats correctly classified as such and the number of beats which triggered a false alarm. Then we tested the four algorithms’ performance in the detection of the three main types of erroneous beats: extra, missed, and misplaced beats. To test the detection of “e” (extra) beats, after inserting a beat before every beat of index *k* = 100 *n* (*n* {1, 2, 3, …}) of each record, we recorded how many of these beats were correctly identified as erroneous beats. Similarly, we measured the performance in the detection of “s” (missed) beats by removing from the original record every beat of index *k* = 100 *n*.

Finally, we assessed each algorithm’s efficacy in the detection of “m” (misplaced) beats, by shifting each beat of index *k* = 100 *n* by a temporal displacement Δ*t*. In order to account for the fact that the detectability of misplaced beats depends on the heart rate variability of the R-R series, we used a record-dependent displacement Δ*t _{r}* =

Table I reports, for the normal series and for each artificially-corrupted series, the number of beats tested (columns “Tot #”) and the percentage of these that were correctly detected (columns labelled “Corr”) by the four algorithms. For the “m” beats, the average displacement $\overline{\mathrm{\Delta}t}={\text{mean}}_{r}\phantom{\rule{0.16667em}{0ex}}\mathrm{\Delta}{t}_{r}$ among all records is also reported.

Results of the detection of irregular beats on artificially-corrupted R-R series. Our method (PP) is compared with three reference methods (Rand 2007, Mateo 2003, and Berntson 1990) and the percentage of correct detection (Corr) of normal (“N”) **...**

As one further goal of our algorithm is to try to correct erroneous beats, it is important that the algorithm is also able to classify the type of error (“e”, “s”, and “m”) and, as a result, the type of action needed to attempt to fix it. All algorithms considered, except the one by Berntson et al., provide this type of information. Table I reports, in the columns labelled “CT”, the percentage of beats for which each algorithm correctly classified the type of error.

These results show that our method represents an improvement over the previous algorithms. In fact it presents the best false detection rate (row “N”) with only 9 normal beats out of 60017 reported as erroneous (which correspond to 99.985% correct classification), while the second best method, by Berntson, triggered a false detection in 37 cases. Our method also shows perfect detection of missed (row “s”) and extra beats (row “e”), and better detection of misplaced beats at all levels of SNR (rows “m” of Table I and Fig. 3). In terms of sensitivity, only the method by Mateo and Laguna provides similar results to our approach but at the cost of a much lower specificity.

The left side of the plot (Corr) reports the percentage of the beats shifted by Δ*t*_{r} = *q* · RMSSD_{r} (as detailed in Section III-A) that were correctly detected as erroneous. The right side of the plot (CT) shows the percentage of the same **...**

We also performed a Monte Carlo analysis of performance by randomly sampling the test beats with probability 1%, obtainining results very similar to those reported here. We decided to report the results of the deterministic sampling of the test beats (every hundredth beat), rather than those from random sampling, to support reproducibility of results and allow other scientists to test and compare their algorithms on exactly the same recordings and beats.

In this section we report the results of evaluating the performance of the correction phase of our algorithm and compare it with existing methods.

We wanted to find out how accurately our model is able to estimate the unknown time of occurrence of a beat (either a missed or a misplaced beat). This was simply evaluated by using (6) to estimate the beat location
${\widehat{u}}_{k}^{\text{PP}}$ of every beat of index *k* = 100 *n* from each original record described in Section III-A and then measuring the difference,
${\widehat{u}}_{k}^{\text{PP}}-{u}_{k}$, between the predicted time of the beat and the actual time of the original beat.

We compared our algorithm with four existing methods of estimating the position of a missing beat. As first method we used the “nguess” function available from the WFDB library of the PhysioToolkit [35] and called its estimate
${\widehat{u}}_{k}^{\text{N}}$. The second method consisted in splitting the interval around the missing beat in two identical R-R intervals by inserting a beat at time
${\widehat{u}}_{k}^{\text{H}}$ such that
${u}_{k+1}-{\widehat{u}}_{k}^{\text{H}}={\widehat{u}}_{k}^{\text{H}}-{u}_{k-1}$. Then we tested the family of * _{N}* estimators given by equation (45) of [36] and found, in agreement with their findings, that the first order estimator (

Table II reports a comparison between the different methods in terms of root-mean-square (RMS) error between the estimated time and the original beat time. The first row shows the result of pooling all 602 beats from all subjects and then evaluating the pooled RMS. The second row reports the result of finding the RMS error for each subject separately and then averaging the results. Similarly, the last row shows the median value of the RMS error for the individual subjects.

Results of the estimation of unknown beat times. The root mean square (RMS) of the error between the beat time estimated by our method
${\widehat{u}}_{k}^{\text{PP}}$ and the true beat time *u*_{k} is compared with that of the other methods reported in the text.

Our algorithm (first column) clearly outperforms all four reference methods in all measures and provides accurate estimates with an RMS error of under 15 ms. This is also confirmed by the probability density function of the difference between the estimated beat time *û _{k}* and the actual time of the original beat

As mentioned before, the algorithm presented here can also be used to detect ectopic events based on the R-R intervals alone, not for the purpose of classifying the different types of arrhythmias but simply to try to replace them with fictitious normal beats and carry on with the desired type of HRV analysis.

The performance of our algorithm in the detection of ectopic beats was evaluated on the MIT-BIH Arrhythmia database [35], [37]. For this analysis we selected the recordings having less than 10 non-normal beats in any 20-second window.^{2} The main reason is that we are interested in R-R series that are mildly affected by erroneous or ectopic events, and that can be recovered and further processed with HRV techniques. In a recording with an excessive number of irregular beats, the information about the underlying sinus rhythm cannot be reliably recovered and therefore the HRV indices cannot be assessed.

For each record of the MIT-BIH Arrhythmia Database considered, we ran the four algorithms described in Section III-A. As some of the algorithms use a 60 s sliding window to adaptively estimate their thresholds, we did not test the detection in the first minute of the record. We also ignored types of ectopic beats that merely affect the morphology of the contraction potential but not its timing because all the algorithm that we tested do not have access to the raw ECG data but only to the series of R-R intervals.

Table III reports the confusion matrices obtained for the four algorithms, whereas Table IV shows their performance measures. Our method outperforms the previous algorithms, especially in terms of a lower number of false detections, resulting in higher specificity and positive predicted value. Figure 5 shows a few examples of real arrhythmic events from this database and how our algorithm detects and corrects them.

The figure shows examples of beats from the MIT-BIH Arrhythmia database and how they were detected and corrected by our algorithm. The markers (little star, triangle, and square) and the associated labels (“N”, “V”, “a”) **...**

Confusion matrices for the point-process method presented here (PP) and for the three state-of-the-art methods considered against the experts’ annotations in the database (Annot). The label “*N*” represent normal beats while “≠ **...**

We present a novel paradigm for detection and correction of ectopic and erroneous heartbeats. At the core of the method is the ability of the underlying probabilistic model to indicate with accurate precision the degree of probability of having a heartbeat at each moment in time. This feature allows to create a tree-like decisional structure based on quantitative comparisons between the observed events and the probability of having the event at that time. More specifically, we test the alternative hypotheses that: a) a beat has been correctly identified, b) a beat has been missed, c) the observed beat is spurious, d) one or two beats need to be moved to a most probable time according to the probability distribution, or e) a beat is a resetting ectopic beat. Since the model is defined in continuous time and is updated recursively in an on-line fashion, the devised algorithm is able to process the decisional tree in real-time using only short-time past information, and then finalizing the correction within at maximum three beats in the future.

Our method, similarly to the IPFM-based method, is based on a physiologically motivated model, with the further advantage that the point process framework allows for dynamically updating the model by tracking the time-varying heartbeat variations, including an objective goodness-of-fit framework which can validate the decisional outcome based on optimal model selection. Importantly, the parameter vector is here computed only to estimate the likelihood, but this is the same vector that defines the instantaneous time and frequency domain point-process HRV as in [17], [38]–[40], thus providing a simultaneous assessment of cardiovascular control that is unbiased by artifactual or arrhythmic events.

Our results provide evidence of the efficacy of the detection method in recognizing missing, spurious, mis-detected, and irregular heartbeats with remarkable accuracy even in scenarios that require discerning very small time resolution errors at millisecond scales. They further point at our method as more accurate than three previously published detection methods and four correction methods. In particular, our automated method achieves 99.985% specificity for detection of correct beats, 100% sensitivity to missed and extra beats, 96.01% sensitivity for beats artificially misplaced by 137 ms on average, and 98.73% positive prediction value for real arrhythmic events. Our algorithm represents an improvement over previous procedures: only one other method provides similar results but at the cost of a much lower specificity. Our correction procedure provides accurate estimates of the replacement time of an unobserved beat, with and error below 15 ms, outperforming other existing methods of estimating the time of missing beats.

One important feature of the technique is that it works with predetermined decisional thresholds that have been fixed once for all during our analyses, and have been demonstrated to work on a wide range of experimental recordings, with high robustness to noise and ability to cope with inter-subject variability. In particular, the algorithm can be used for exclusive pre-processing purposes of any R-R interval series and consequent application of any kind of HRV analysis. The routine is very fast, easy to run, and available online at no cost for research purposes [19]. We hope dissemination will further corroborate the robustness and efficacy of the proposed method in an even wider range of physiological and clinical settings.

The authors would like to thank Prof. Roger Mark for his valuable comments and suggestions that were of great help during the first stages of this work.

This work was supported by National Institutes of Health (NIH) through grants R01-HL084502 (R.B.) and DP1-OD003646 (E.N.B.).

**Luca Citi** (S’06–M’09) was born in Pistoia, Italy, in 1979. He received a laurea (MS) degree in Electronic Engineering with major in Biomedical Engineering, in 2004 from Università degli Studi di Firenze (Florence, Italy). In 2009, he received his PhD in Biorobotics Science and Engineering jointly offered by IMT Lucca (Lucca, Italy) and Scuola Superiore Sant’Anna (Pisa, Italy).

He is currently a post-doctoral research fellow at Massachusetts General Hospital/Harvard Medical School and affiliated with the Massachusetts Institute of Technology (Cambridge, MA, U.S.A.) where he is part of the “Neuroscience Statistics Research Laboratory” led by Prof. Emery Brown, and the “Neuro–Cardiovascular Signal Processing Unit” led by Dr. Riccardo Barbieri. He works on the statistical analysis of human cardiovascular rhythms and of neural spike trains. His work on HRV focuses on devising and refining stochastic models of cardiovascular control where heart beat intervals are governed by a history-dependent point process. Other interests include signal processing in computational neuroscience and neural interfaces.

**Emery N Brown** (M’01, SM’06, F’08) received the B.A. degree from Harvard College, Cambridge, MA, the M.D. degree from Harvard Medical School, Boston, MA, and the A.M. and Ph.D. degrees in statistics from Harvard University, Cambridge, MA. He is currently a Professor of computational neuroscience and health sciences and technology in the Department of Brain and Cognitive Sciences and the Division of Health Sciences and Technology, Harvard-Massachusetts Institute of Technology, the Warren M. Zapol Professor of Anesthesia at Harvard Medical School, and an Anesthesiologist in the Department of Anesthesia, Critical Care and Pain Medicine, Massachusetts General Hospital. His research interests include signal processing algorithms for the study of neural systems, functional neuroimaging, and electrophysiological approaches to study general anesthesia in humans. Dr. Brown is a member of the Association of University of Anesthesiologists. He is also a Fellow of the American Statistical Association, a Fellow of the American Association for the Advancement of Science, a member of the Institute of Medicine, and a Fellow of the American Academy of Arts and Sciences. He was a recipient of a 2007 National Institutes of Health Director’s Pioneer Award and the 2011 recipient of the National Institute of Statistical Sciences Sacks Award for Outstanding Cross-Disciplinary Research.

**Riccardo Barbieri** (M’00, SM’08) was born in Rome, Italy, in 1967. He received the M.S. degree in Electrical Engineering from the University of Rome “La Sapienza”, Rome, Italy, in 1992, and the Ph.D. in Biomedical Engineering from Boston University, Boston, MA, in 1998.

Dr. Barbieri is currently Assistant Professor of Anaesthesia at Harvard Medical School, Massachusetts General Hospital, Research Affiliate at Massachusetts Institute of Technology and Principal Investigator of the Neuro-Cardiovascular Signal Processing Unit within the Neuroscience Statistics Research Laboratory. His main research interests include: modeling of biological systems (with emphasis on the cardiovascular system, the central nervous system and the autonomic nervous system); biomedical signal processing (with emphasis on cardiovascular control and cardiovascular variability signals); and signal processing in computational neuroscience (EEG, Evoked Potentials, single cell neuronal recordings).

In this Appendix we find the probability distribution of *τ* − *u _{k}* which is the sum of two unknown intervals

$$\tau -{u}_{k}={\widehat{w}}_{k+1}+{\widehat{w}}_{k+2}={\mu}_{1}+{\epsilon}_{1}+{\mu}_{2}({\mu}_{1})+{\stackrel{\sim}{\theta}}_{1}({u}_{k}){\epsilon}_{1}+{\epsilon}_{2}$$

(7)

which leads to *μ*_{1+2} = E[*τ* − *u _{k}*] =

Following a similar reasoning, we obtain that the probability distribution of the sum of three unknown intervals can be modelled as an IG distribution with parameters *μ*_{1+2+3} = *μ*_{1} + *μ*_{2}(*μ*_{1}) + *μ*_{3}(*μ*_{1}, *μ*_{2}) and

$${\lambda}_{1+2+3}=\frac{{\lambda}_{1}{\mu}_{1+2+3}^{3}}{{(1+{\stackrel{\sim}{\theta}}_{1}+{\stackrel{\sim}{\theta}}_{2})}^{2}{\mu}_{1}^{3}+{(1+{\stackrel{\sim}{\theta}}_{1})}^{2}{\mu}_{2}^{3}({\mu}_{1})+{\mu}_{3}^{3}({\mu}_{1},{\mu}_{2})}.$$

^{1}The records used were: f1o01, f1o05, f1o10, f1y01, f1y03, f1y08, f2o04, f2y03, f2y06.

^{2}The records used were: 100, 101, 103, 105, 108, 112, 113, 114, 115, 116, 117, 121, 122, 123, 215, 230.

1. Camm A, Malik M, Bigger J, Breithardt G, Cerutti S, Cohen R, Coumel P. Heart rate variability: standards of measurement, physiological interpretation and clinical use. task force of the european society of cardiology and the north american society of pacing and electrophysiology. Circulation. 1996;93(5):1043. [PubMed]

2. Acharya U, Joseph P, Kannathal N, Lim C, Suri J. Heart rate variability: a review. Med Biol Eng Comput. 2006;44:1031–1051. [PubMed]

3. Acharya U, Suri J, Spaan J, Krishnan S. Advances in cardiac signal processing. Springer Verlag; 2007.

4. Lippman N, Stein KM, Lerman BB. Comparison of methods for removal of ectopy in measurement of heart rate variability. Am J Physiol — Heart and Circulatory Physiology. 1994 Jul;267(1):H411–H418. [PubMed]

5. Kamath MV, Fallen EL. Correction of the heart rate variability signal for ectopics and missing beats. In: Malik M, Camm AJ, editors. Heart rate variability. 1995. pp. 75–85.

6. Clifford G, McSharry P, Tarassenko L. Characterizing artefact in the normal human 24-hour RR time series to aid identification and artificial replication of circadian variations in human beat to beat heart rate using a simple threshold. Proceedings of Computers in Cardiology; Memphis. Sep. 2002; pp. 129–132.

7. Cheung MN. Detection of and recovery from errors in cardiac interbeat intervals. Psychophysiology. 1981 May;18(3):341–346. [PubMed]

8. Malik M, Cripps T, Farrell T, Camm A. Prognostic value of heart rate variability after myocardial infarction. a comparison of different data-processing methods. Med Biol Eng Comput. 1989;27(6):603–611.

9. Berntson GG, Quigley KS, Jang JF, Boysen ST. An approach to artifact identification: Application to heart period data. Psychophysiology. 1990 Sep;27(5):586–598. [PubMed]

10. Weber EJM, Molenaar PCM, Molen MW. A nonstationarity test for the spectral analysis of physiological time series with an application to respiratory sinus arrhythmia. Psychophysiology. 1992 Jan;29(1):55–62. [PubMed]

11. Sapoznikov D, Luria M, Mahler Y, Gotsman M. Computer processing of artifact and arrhythmias in heart rate variability analysis. Computer methods and programs in biomedicine. 1992;39(1–2):75–84. [PubMed]

12. Brennan M, Palaniswami M, Kamen P. A new model-based ectopic beat correction algorithm for heart rate variability. Proceedings of the 23rd Annual International Conference of the IEEE EMBS; 2001. pp. 567–570.

13. Mateo J, Laguna P. Analysis of heart rate variability in the presence of ectopic beats using the heart timing signal. IEEE Trans Biomed Eng. (3) 2003 Mar;50:334–343. [PubMed]

14. Barbieri R, Brown EN. Correction of erroneous and ectopic beats using a point process adaptive algorithm. Proceedings of the 28th Annual International Conference of the IEEE EMBS; Sep. 2006; pp. 3373–3376. [PubMed]

15. Rand J, Hoover A, Fishel S, Moss J, Pappas J, Muth E. Real-Time correction of heart interbeat intervals. IEEE Trans Biomed Eng. 2007 May;54(5):946–950. [PubMed]

16. Widjaja D, Vandeput S, Taelman J, Braeken M, Otte R, Van den Bergh B, Van Huffel S. Accurate R peak detection and advanced preprocessing of normal ECG for heart rate variability analysis. Proceedings of Computing in Cardiology; Belfast. Sep. 2010; pp. 533–536.

17. Barbieri R, Matten EC, Alabi AR, Brown EN. A point-process model of human heartbeat intervals: new definitions of heart rate and heart rate variability. Am J Physiol — Heart and Circulatory Physiology. 2005;288:H424–435. [PubMed]

18. Citi L, Brown EN, Barbieri R. A point process local likelihood algorithm for robust and automated heart beat detection and correction. Proceedings of Computing in Cardiology; Hangzhou. Sep. 2011.

19. Citi L, Barbieri R. Point process detection and correction of erroneous and ectopic heart beats. 2012 [Online]. Available: http://users.neurostat.mit.edu/barbieri/ncspu. [PMC free article] [PubMed]

20. Moody G. Spectral analysis of heart rate without resampling. Proceedings of Computers in Cardiology, IEEE; 1993. pp. 715–718.

21. Guyton A, Hall J. Textbook of medical physiology. Vol. 9 Saunders Philadelphia, PA: 1991.

22. Hyndman B, Mohn R. A model of the cardiac pacemaker and its use in decoding the information content of cardiac intervals. Automedica. 1975;1:239–252.

23. DeBoer R, Karemaker J, Strackee J. Description of heart-rate variability data in accordance with a physiological model for the genesis of heartbeats. Psychophysiology. 1985;22(2):147–155. [PubMed]

24. Berger R, Akselrod S, Gordon D, Cohen R. An efficient algorithm for spectral analysis of heart rate variability. IEEE Trans Biomed Eng. 1986;9:900–904. [PubMed]

25. Mateo J, Laguna P. Improved heart rate variability signal analysis from the beat occurrence times according to the ipfm model. IEEE Trans Biomed Eng. 2000;47(8):985–996. [PubMed]

26. Brown E. Methods and Models in Neurophysics. Vol. 14. Elsevier; 2005. Theory of point processes for neural systems, ser; pp. 691–726.

27. Feng J. Computational neuroscience: a comprehensive approach. CRC press; 2004.

28. Chhikara R, Folks L. The inverse Gaussian distribution: theory, methodology, and applications. Vol. 95 CRC; 1989.

29. Gerstein G, Mandelbrot B. Random walk models for the spike activity of a single neuron. Biophysical Journal. 1964;4(1):41–68. [PubMed]

30. Tuckwell H. Introduction to theoretical neurobiology: Nonlinear and stochastic theories. Vol. 2 Cambridge Univ Pr; 2005.

31. Stanley G, Poolla K, Siegel R. Threshold modeling of autonomic control of heart rate variability. IEEE Trans Biomed Eng. 2000;47(9):1147–1153. [PubMed]

32. Andersen P. Statistical models based on counting processes. Springer Verlag; 1993.

33. Loader C. Local regression and likelihood. Springer Verlag; 1999.

34. Iyengar N, Peng C, Morin R, Goldberger A, Lipsitz L. Age-related alterations in the fractal scaling of cardiac interbeat interval dynamics. Am J Physiol — Regulatory, Integrative and Comparative Physiology. 1996;271(4):R1078–R1084. [PubMed]

35. Goldberger A, Amaral L, Glass L, Hausdorff J, Ivanov P, Mark R, Mietus J, Moody G, Peng C, Stanley H. Physiobank, physiotoolkit, and physionet: Components of a new research resource for complex physiologic signals. Circulation. 2000;101(23):e215–e220. [PubMed]

36. Solem K, Laguna P, Sörnmo L. An efficient method for handling ectopic beats using the heart timing signal. IEEE Trans Biomed Eng. 2006;53(1):13–20. [PubMed]

37. Moody GB, Mark RG. The impact of the MIT-BIH arrhythmia database. IEEE Eng Med Biol Mag. 2001 Jun;20(3):45–50. [PubMed]

38. Barbieri R, Brown E. Analysis of heartbeat dynamics by point process adaptive filtering. IEEE Trans Biomed Eng. (1) 2006;53:4–12. [PubMed]

39. Chen Z, Brown E, Barbieri R. Assessment of autonomic control and respiratory sinus arrhythmia using point process models of human heart beat dynamics. IEEE Trans Biomed Eng. 2009;56(7):1791–1802. [PMC free article] [PubMed]

40. Chen Z, Brown E, Barbieri R. Characterizing nonlinear heartbeat dynamics within a point process framework. IEEE Trans Biomed Eng. 2010;57(6):1335–1347. [PMC free article] [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |