Method and Apparatus for Detecting Paroxysmal Atrial Fibrillation
This invention relates to a method and apparatus for detecting probable paroxysmal atrial fibrillation.
Atrial fibrillation (AF) is the most common heart arrhythmia in clinical practice, and has serious morbidity and mortality. AF is a significant risk factor for stroke; about 15% of strokes occur in people with AF. The prevalence of AF increases with age, with an estimated 10 % of the 80-89 age group suffering from AF. AF can either be chronic or intermittent; intermittent AF is referred to as paroxysmal atrial fibrillation (PAF).
Atrial fibrillation (AF) is currently diagnosed by searching for episodes of atrial fibrillation detected using a 12 lead electrocardiogram, an ambulatory Holter monitor (which records a single or two-channel ECG over 24-72 hours) or an event recorder (which records one or two channels of ECG when triggered by patient). An episode of AF can be recognised by the presence of a low-amplitude variable signal prior to ventricular depolarisation and the absence of normal amplitude P waves. If an episode of AF occurs during the recording time the patient is diagnosed as suffering from atrial fibrillation, otherwise not. No attempt is currently made to diagnose atrial fibrillation in cases where no actual episode of AF is recorded.
US2002143266 discloses a technique for detecting episodes of atrial fibrillation based on analysis of R- intervals and beat types. The technique uses a hidden Markov model as a classification technique, and also carries out detection of P waves to assist in classification. Similarly, WO0224068 describes a technique for using the R-R difference series to differentiate between sections of electrocardiogram where a subject is in normal rhythm and in atrial fibrillation. US6178347 discusses another generic episode-based classification system which subtracts out the QRS complex, and then carries out spectral analysis of the residual signal in order to identify episodes of atrial fibrillation and atrial flutter. This idea is also known in the technical literature (for example, see Stridh M. and L. Sornmo. Spatiotemporal QRST cancellation techniques for analysis of atrial fibrillation. IEEE Trans on Biomed Eng, 48: 105-111, January 2001). However all the above work is based on detecting episodes of atrial fibrillation, and makes no attempt to assess the likelihood of atrial fibrillation using the non-AF-episode information.
However, some prior art exists in the assessment of RR interval dynamics in the immediate period prior to the onset of an episode of atrial fibrillation. Nikman et al. (Nikman S., T Makikallio, S. Yli-Mayry, S. Pikkujamsa, A.M. Koivisto, P. Reinikainen, K.E. Johani Airaksinen, H.N. Huikuri. Altered complexity and correlation properties of RR interval dynamics before the spontaneous onset of paroxysmal atrial fibrillation. Circulation; 100: 2079-2084, 1999) noted a change in the R-R dynamics in pre-episode ECG as assessed by a measure called approximate entropy. The present invention differs by using spectral features and detection of atrial premature contractions (rather than approximate entropy) to derive information from the non-AF-episode ECG. Previous work has also considered assessment of autonomic tone using heart rate variability to assess the likelihood of an AF episode being initiated (Huang J.L., Z.C. Wen, .L. Lee, M.S. Chang, S.A. Chen. Changes in autonomic tone before the onset of paroxysmal atrial fibrillation. International J of Cardiology, 66:275-283, 1998., Hickey B., C. Heneghan. Screening for paroxysmal atrial fibrillation using atrial premature contractions and spectral measures. Computers in Cardiology, 29:217-220, 2002. De Chazal P., Heneghan, C. Automated assessment of atrial fibrillation. Computers in Cardiology, 28:117-120, 2001). However, the present invention differs in that it uses spectral measures to initially divide the classification problem into two subclasses, and it does not solely focus on predicting AF onset. Finally,, some work exists in the area of using atrial premature contractions as a possible marker for the onset of AF (Waktare J.E.P., K. Hnatkova, S.M. Sopher, F.D. Murgatroyd, X. Guo, A J. Cani , M. Malik. The role of atrial ectopics in initiating paroxysmal atrial fibrillation. European Heart Journal, 22:333-339, 2001, Hickey B., C. Heneghan. Screening for paroxysmal atrial fibrillation using atrial premature contractions and spectral measures. Computers in Cardiology, 29:217-220, 2002). However, the present invention differs in that it combines information from atrial premature contractions together with spectral R-R measures, and also does not solely focus on predicting AF onset. It is an object of the invention to provide a method and apparatus for detecting probable paroxysmal atrial fibrillation (PAF) in a subject from a prolonged recording of an electrocardiogram (ECG) of that subject not necessarily containing an episode of AF.
Accordingly, the invention provides a method of detecting probable paroxysmal atrial fibrillation (PAF) in a subject from a prolonged recording of an electrocardiogram (ECG) of that subject not necessarily containing an episode of atrial fibrillation, the method comprising the following steps: (A) for each of a plurality of segments of the ECG: (i) determining cardiac intervals, (ii) calculating the correlation between the high and low frequency spectral content of the cardiac interval sequence at predetermined intervals during the ECG segment, and (iii) classifying the subject into a first or second group according to whether the degree of correlation is above or below a predetermined threshold, if the subject is in the first group: (iv) determining the power spectral density (PSD) of a cardiac interval sequence at a plurality of frequencies for at least a part of the segment, and (v) classifying the magnitudes of the power at each frequency to provide a probability of PAF, if t s sυbjed is in 1be second group: (vi) determining the number of atrial premature contractions (APCs) in the segment, and, if the number of APCs in the segment exceeds a predetermined amount, (vii) determining the power spectral density (PSD) of a cardiac interval sequence at a plurality of frequencies for at least a part of the segment, and (viii) classifying the magnitudes of the power at each frequency to provide a probability of PAF, the method further including: (B) providing a positive indication of PAF in a subject as a function of the PAF probabilities provided by steps (v) and (viii) for all of the said plurality of ECG segments. The invention further provides an apparatus adapted to perform the above method and a computer-readable storage medium containing instructions which, when executed by a computer, will perform the above method.
In the present specification the term "cardiac interval" means a cardiac interbeat interval, such as the RR, PP, TT interval or other time interval between electrically equivalent points on the cardiac cycle, or a cardiac intrabeat interval such as the PR or QT interval. Although the embodiment to be described uses the same RR cardiac interbeat interval for each of steps A(ii), A(iv) and A(vii), this is not essential and the same or a different cardiac interval could be used for each of these steps.
This invention ignores the conventional approach of diagnosing paroxysmal atrial fibrillation. Conventional methods depend heavily on whether episodes of AF occur during the ECG recording time, whereas the current invention does not. If no episode occurs during the recording time, the test fails using conventional methods. However, with this invention, the test can still succeed. Furthermore, this invention can be based on a single lead ECG, although not limited thereto, unlike most current techniques that use multiple leads.
An embodiment of the invention will now be described,, by way of example only, with reference to the single figure of the drawings which is a flow diagram of the preferred embodiment.
Step 1 : ECG Acquisition
An electrocardiogram (ECG) of the subject was recorded over a 10 hour period. This was done using a 12 lead ECG system. For the present embodiment only a single lead needs to be used. This may be, for example, lead II or a modified lead V5. Alternatively an ambulatory Holter monitor or an event recorder may be used to obtain the ECG.
Step 2: ECG Pre-processing
The raw ECG signal acquired as above was processed using a linear phase band pass FIR filter windowed using a Kaiser-Bessel window to reduce baseline wander and EMG artefact. The band stop frequencies were set at 0.25 Hz and 20 Hz. QRS detection was implemented using a Hubert Transform method. From this set of R peaks, the RR intervals were determined by known techniques.
The ECG was now divided up into twenty contiguous 30-minute segments and the following two-stage classification system performed for each segment.
Steps 3 - 5: Classification Stage 1
Since changes in autonomic tone have been implicated in the initiation of AF, the first stage of classification provides a measure of such changes. In this first stage, the correlation between the high frequency (HF) and low frequency (LF) spectral content of the RR interval sequence was calculated at 2-minute intervals during the 30 minute segment. These components were calculated as follows.
Step 3 : The HF and LF components were derived from the interval based power spectral density (PSD) estimates of non-overlapping two-minute RR interval sequences. The RR interval sequence contains both normal and ectopic beats, i.e., no attempt was made to restrict the analysis to NN intervals. To obtain a zero-mean sequence the mean RR interval was subtracted from each segment. The sequence was zero padded to the nearest integer multiple of 98 exceeding the length of the sequence, and the Fast Fourier Transform (FFT) was taken of the entire sequence. The absolute values of the FFT coefficients were then squared to yield a periodogram estimate of the PSD. Adjacent frequency bins were then combined to result in a 98 point PSD estimate (of which only bins 0-49 are relevant since 50-97 provide identical information as 1-48). The magnitudes of these PSD bins from 0.04 to 0J5 cycles/interval were added to yield the LF component, and similarly the bins from 0J5 to 0.4 cycles/interval yielded the HF component (this is in line with the conventional frequency ranges used in defining LF and HF components).
Step 4: For the 30-minute segment, this process yielded two sets of 15 values, and the correlation coefficient between these sets of numbers was taken, where the correlation coefficient was calculated as follows:
where x
( and y
t are the LF and HF powers respectively, x and y are their mean values, and σ
x and σ
y are the standard deviations.
Step 5: From this analysis, two subgroups of records were developed - those records with a correlation coefficient above 0.75 (group A), and those records with a correlation coefficient below 0.75 (group B). A physiological interpretation is that the subjects in Group A experience parallel changes in sympathetic and parasympathetic activation, whereas those in Group B experience some decoupling. This decoupling could be due to something as simple as a postural change, or may reflect an origin of unknown cause.
Steps 6 - 10: Classification stage 2
As a result of stage 1, there are now two groups^ Group A and Group B, both of which contain both PAF and non-PAF patients. To separate the PAF patients from the non-PAF patients in each group, further heart rate variability (HEN) analysis was carried out.
To separate the members of Group A from each other, the PSD of the RR interval sequence was once again taken and calculated as above. Adjacent frequency bins were then combined to result in a 64 point PSD estimate, Step 7, of which only bins 0-32 are relevant since 33-63 provide identical information as 1-31. This analysis was carried out on the final 10 minutes of the 30 minute segment. The magnitudes of these PSD bins were used as features in a linear discriminant classifier (LDC), Step 7, to be described below.
To separate the members of group B from each other, first the number of Atrial Premature Contractions (APCs) in the 30-minute segment was determined, Step 8, using known techniques. Such determination may also include analysing changes in the shape or magnitude of the P wave, QRS complex or T wave. If there were zero or one APC in the segment, the probability PPAF (0<PPAF<1) of the subject having PAF is very small, and is here arbitrarily given a figure of 0.01.
If the APC count was greater than one, the same analysis was carried out in Step 9 as was carried out in Step 6 for Group A subjects, but this time over the final 20 minutes of the 30-minute segment. Again, Step 10, the magnitudes of the PSD bins were used as features in a linear discriminant classifier.
For Steps 7 and 10 a linear discriminant classifier (LDC), based on Fisher's rule, was used. For each 30-minute segment of data fed in, there are two possible output classes - PAF or non-PAF. The classifier outputs a set of numbers representing the probability estimate of each class, in response to a set of input features. Linear discriminants partition the feature space into different classes using a set of hyper- planes. Optimisation of the model is achieved through direct calculation and is extremely fast relative to other models such as neural networks. The training of a LDC proceeds as follows. Let s be a dx 1 column vector containing feature values calculated from a data set (e.g., the mean RR interval, the standard deviation of the RR intervals, etc.). We wish to assign ∑ to one of c possible classes (c=2 in our case). A total of N feature vectors are available for training the classifier, with the number of feature vectors representing class k equal to N/c., i.e.:
N = ∑Nlc (2) k The nth training vector in class k is denoted as __k . The class-conditional mean vectors γ_k are defined as:
We now define a common covariance matrix defined over all classes (i.e., we assume that each class only differs in its mean value, and not in its higher order statistics). The common covariance matrix is defined as:
∑ = - Nτ —-c ∑k=l _n£=ϊ (**„~ Vk k,n - μ* ϊ (4)
Given these definitions, we can now calculate a discriminant valuer for an arbitrary data vector : yk = -iμ*∑~ fc + μ*∑_lχ + losO/. ) (5) where u is the a priori probability of the vector x being from class k. It is easy to convert the discriminant values to posterior probabilities using:v{ π) _ ______ (6) ∑exp(yk)
This formulation provides a mapping from discriminant value to posterior probabilities wliich tends to bias the probabilities towards the extremes (since linear discriminant classifiers generally tend to underestimate the class probabilities). The final class assigned to x is the class with the highest posterior probability.
Two separate LDC models (with differing parameters) were used: one to separate the PAF and non-PAF members of Group A at Step 7, and a second to separate the PAF and non-PAF members of Group B at Step 10. Sensitivity and specificity can be traded off by adding terms to Eq. (5) which weight the cost of assigning data to the wrong class, or equivalently by altering the threshold at which we make a decision to assign a vector to a class based on posterior probability. For the present embodiment, the posterior probability of the PAF class had to exceed 0.6 (60%) for a 30-minute segment to be classed as PAF in origin.
For both LDC models, an exhaustive search was performed to find the best two spectral features that optimally separated PAF members from non-PAF members. To reduce computational complexity, only the best 2-feature combination was calculated. For Group A, the two optimal spectral measures corresponded approximately to spectral features at frequencies of 0J 1 cycles/interval and 0.48 cycles/interval, whereas for Group B the two optimal spectral measures corresponded approximately to spectral features at frequencies of 0.016 cycles/interval and 0.375 cycles/interval.
The above analysis outputs a probability PPAF for each 30-minute segment. A 30- minute recording was classified as PAF in origin if PPAF was greater than or equal to 60%. However, in order to make a decision on 10 hours of data, the ECG was classified as coming from a AF patient if:
(a) Over 10 % of the segments were classified as PAF (PPAF >= 60%), and (b) the average probability of PAF ( PPAF) over all segments was greater than 35%.
The above method was tested on ECG records drawn from three databases:
The Physionet Normal Sinus Rhythm Database (NSRD). The database consists of 17 records, each approximately 24 hours long, recorded using lead II. These patients have no significant arrhythmias, and include 5 men, aged 26 to 45, and 13 women aged 20 to 50.
Measurements taken at University College Dublin (UCDDB) using a Reynolds ambulatory Holter monitor. The database consists often overnight records of approximately 8-hours duration from subjects with no known cardiac pathology, aged between 23 and 40 (modified N5). The signals were sampled at 128 Hz, with 12-bit resolution.
The MIT-BIH AF database (AFDB); This database includes 23 long-term ECG recordings of human subjects with atrial fibrillation (mostly paroxysmal). The individual recordings are each approximately 10 hours in duration, and contain two ECG signals each sampled at 250 Hz with 12-bit resolution. The number of episodes of AF varies from 0 to 39.
The method was tested on the ΝSRDB and UCDDB, in order to check the specificity, and then on the AFDB, in order to test the sensitivity. The results are shown in Tables 1 to 3.
ΝSRDB
*at least 10% of the recording is greater than 60% PAF
Table 1: Normal Sinus Rhythm Database
UCDDB
*at least 10% of the recording is greater than 60% PAF
Table 2: University College Dublin Database AFDB
*at least 10% of the recording is greater than 60% PAF
Table 3: Atrial fibrillation database
The classification system achieved 100% accuracy in separating 26 Non-PAF and 23 PAF patients.
Apart from the recording of the initial ECG, Step 1, it will be recognised by the person skilled in the art that the entire method can be implemented in software on a programmable general computer, such as a PC. However, this does not rule out the method being implemented in special purpose devices in, for example, microcode. Further, while the ECG used in the above embodiment is taken externally of the subject, it is quite possible to obtain the ECG from implanted devices.
Variations of the above method are possible. For example, it is not necessary to use 10 hours of ECG recording, but preferably at least 5 hours should be used. Further, the individual segments subject to Steps 2 to 10 need not be 30 minutes ling, but they should preferably be at least 15 minutes long. Also, the segments need not be contiguous — they could overlap (slower processing but greater accuracy) or be non-contiguous (faster processing but reduced accuracy). It is also not necessary that the RR intervals be 2 minutes in Step 3, but they are preferably at least 1 minute intervals. Furthermore, the portion of each 30-minute segment used in Steps 6 and 9 could be the whole segment, but preferably in each case it is a significant portion of the segment at or near the end of the segment. Also, the number of frequencies used in Steps 6 and 9 need not be 64 but is preferably at least 16. Also, in Step 2, although the RR intervals were determined for use in Step 3, a different cardiac interval could be used, for example, the RR, PP, TT interval or other time interval between electrically equivalent points on the cardiac cycle. Alternatively, a cardiac intrabeat interval such as the PR or QT interval could be used. Furthermore, although the embodiment used the RR cardiac interbeat interval for each of Steps 3, 6 and 9, this is not essential and the same or a different cardiac interbeat or intrabeat interval could be used for each of these steps.
The invention is not limited to the embodiment described herein which may be modified or varied without departing from the scope of the invention.