Ventricular tachycardia (VT) occurs most commonly in patients with structural heart disease and can be associated with an increased risk of sudden cardiac death (SCD)  . Conventional techniques such as Implantable Cardioverter-Defibrillators (ICDs) are effective treatments for cardiomyopathy, and have been shown to reduce cardiovascular death and total mortality in patients with reduced ejection fraction after myocardial infarction  . However, ICDs play a limited role in that they only terminate VTs, not prevent them. Furthermore, ICD shocks can be painful for the patient; they increase mortality, and worsen quality of life  . Catheter ablation therapy has progressed significantly over the past few years and is now an important option to control recurrent VTs  . Catheter ablation can prevent recurrent VT without the need for long-term antiarrhythmic therapy, and reduce the need for ICD shocks  .
Current ablation techniques focus on isolating scar tissue by targeting reentry circuit isthmuses and their exits  . To accomplish this, the region of scar must first be defined with the assistance of an electroanatomical system. Substrate mapping can be performed using voltage/pace mapping, followed by registration of electrogram features suggestive of abnormal myocardium. Though passive substrate mapping alone can be used to guide the ablation therapy, combining it with short episodes of activation, and entrainment mapping to isolate slow conduction channels within myocardium increases the likelihood of successful ablation  . However, such methods have a limited success rate, with approximately 50% recurrence rate within 2 years  .
Programmed stimulation for investigation and induction of arrhythmia was initially applied to atrioventricular re-entrant arrhythmia   . The same technique was subsequently applied to the ischemic substrate in the ventricle  and formed the basis of programmed electrical stimulation for ventricular arrhythmia. It has long been recognized that close-coupled pacing, which is achieved through programmed stimulation, can alter local tissue electrophysiologic properties, thus resulting in a change in action potential duration, conduction velocity and refractory periods  . Close-coupled pacing around a scar has been found to induce local electrogram fractionation. This phenomenon is discrete, and observed only at areas of abnormal tissue  . This local slowing of conduction with paced extrastimuli is the electrophysiologic basis of re-entry and as such represents a risk of arrhythmia.
Fractionation increases as tissue is stressed by increased rate or premature beats  . Programmed stimulation is designed to create local tissue conduction block (tissue reaches effective refractory period and therefore cannot propagate signal) in susceptible tissue in order to induce arrhythmia   . The paced electrogram fractionationation analysis (PEFA) protocol was first introduced by Saumarez et al. and has been particularly enlightening in demonstrating an increase in fractionation with prematurity of the signal. When many electrodes are used, the presence of fractionation correlates very well with arrhythmia risk in a prospective fashion  . Limiting the effectiveness of this approach is the time consuming nature of PEFA, increasing associated procedural risks and negatively impacting the willingness of electrophysiologists to perform such a protocol. Therefore, there is a lack of clinical evidence to support the hypothesis of benefit in targeting of ablation. In this paper we propose an automated version of a modified PEFA protocol for VT patients, which expands on previously defined fractionation features by Saumarez et al., and which substantially decreases associated study times.
2. Materials and Methods
Data were obtained from six patients (age 68 ± 12 years, 4 male). All patients were diagnosed with either hypertrophic cardiomyopathy or an arrhythmia (or both in the case of two patients) in the ventricles. Intracardiac bipolar electrogram segments were recorded (sampling rate 1200 Hz) using a 3-dimensionalnavigation system (EnSite™ NavX, St. Jude Medical, MN).
A diagnostic catheter (Inquiry™, St. Jude Medical) was positioned at the right ventricular apex (RVA) as our pacing catheter. A multipolar ablation catheter (Therapy™ Cool Path Duo™, St. Jude Medical) was positioned at various locations in the left ventricle (LV) of two patients, and the right ventricle (RV) of four patients as our roving catheter. An example of our catheter placement during data collection in the RV can be seen in Figure 1 (bottom). Signals were recorded during programed stimulation.
The effective refractory period of the RVA (VERP) was determined using a standard pacing protocol with 10 ms decrement. Thereafter the experimental protocol was performed prior to ablation during native rhythm from the RVA using a pacing sequence which consisted of a paced train (S1) at basic cycle length of 600 ms and an extrastimulus (S2) applied at every seventh beat. The S2 coupling intervals were normalized by the VERP for each patient and three S2 data sets were recorded at each location. These consisted of VERP +150 ms, VERP +100 ms and VERP +50 ms. Figure 2 illustrates the pacing protocol. These data sets were stored for subsequent offline analysis described below. Paced train cycle length was reduced to 500 ms in the presence of frequent spontaneous ventricular ectopy.
Fractionation in our study was defined by conduction time (CT), electrogram duration (ED), and number of deflections (ND). CT refers to the duration between stimulus and onset of the local activation. ED was defined as the duration of the local bipolar electrogram activation, with ND being a measure of the number of deflections (extrema) measured within this activation. An example of the features extracted from the electrogram can be found in Figure 3.
An operator manually annotated these features using NavX system. Annotation was done using the distal channel (D-2) for both the roving and the pacing electrograms. Figure 1 (top) shows a screen shot of the system used and the methodology involved in annotation. Our process of automatically isolating fractionation can be seen in the flowchart in Figure 4. The following section will describe in more detail the steps outlined in this figure. All data analysis was done in Matlab (Mathworks Inc).
Preprocessing consisted of a 4th order Butterworth low pass filter (with a 3 dB cutoff
Figure 1. (Top) An example the interface of the NavX system where manual annotation was done. Here, the MAP channels represents the electrograms from the roving catheter while the RV electrogram represents the electrograms from the RV pacing catheter. Annotation was done using the distal channel (D-2) for both the roving and the pacing electrograms. Here, the caliper (vertical red lines) is expanded to show the electrogram duration (of 56 ms) for the response spike for the fourth stimuli in the drive train. (Bottom) Artist’s rendition of catheter placement in the heart during RV data collection. The RV pacing catheter was used for pacing and the roving catheter was used for data collection.
Figure 2. A processed signal from the RV pacing catheter showing three automatically annotated pacing sequences. With paced train (S1) consisting of 6 beats our extrastimuli (S2) consisting of 1 beat.
Figure 3. Signal from RV channel (green line) is superimposed onto signal from Roving channel (blue line). Three features extracted from electrogram: Conduction time (CT), EGM duration (ED), and number of deflections (ND), are shown. Onset of CT begins at identified pace spike. All pertinent peaks have been numbered for a total ND of 4.
Figure 4. Flowchart showing the progression of our algorithm. The decisional block labeled “Interference” identifies possible interference due to sinus or ectopic beats while the decisional block labeled ‘Physiologically Sound’ checks if isolated CT fits physiological constraints using relative positions of the RV catheter and roving catheter to determine velocity.
frequency of 300 Hz) to remove high frequency disturbances and Teagers-Kaiser Energy operator to simplify peak detection (Figure 5). Teager-Kaiser Energy operator (TK) is a nonlinear quadratic operator proposed by Teager  to measure the real physical energy of a system (En):
where xn is a discrete time signal and n is the discrete time index. This operator is capable of tracking instantaneously-varying special patterns making it ideal for isolating and emphasizing peaks in intracardiac recordings. Finally, signals with low variance are discarded to account for insufficient catheter contact.
2.2. Isolating Pacing Peaks
Signals obtained from the RV pacing catheter are normalized between 0 and 1 to
Figure 5. Original (solid) and processed with Teager-Kaiser energy operator (dotted).
faciletate thresholding of peaks of amplitude less than 0.2; a value which was determined through expert opinion and testing. Furthermore, all peaks which occurred within 250 ms (VERP is rarely less than 200 ms  ) of their neighbor were removed. Our algorithm separates the pacing train and extra stimulibased on the stipulation that the S2 cycle length is always shorter relative to S1; and the fact that the paced train always consists of a minimal of two stimuli. Figure 2 shows an example of automatically isolated pacing peaks from the RV pacing channel. This is an example of a clean signal where no auxiliary cardiac activity had to be accounted for. Unfortunately, this is not always the case. Ectopic beats consist of an irregular beat caused by activation with an origin outside of the cardiac conduction system while sinus beats consist of the heart’s natural contractions. Such physiologically induced interference can make the pacing stimulus inert by previously depolarizing the ventricle. Or it can induce error in our feature metrics if an isolated activation is attributed to our stimulus as opposed to the interference. To account for this, the surface ECG is used.
2.3. Matched Filtering and Non-Pacing Interference
Activations in the surface ECG are isolated using matched filtering  . The template we use consists of a single pacing activation identified on the surface ECG by an expert. Ectopy and fused beats do not commonly follow standard activation, and therefore predominately have no corresponding activation in the surface ECG using matched filtering. The corresponding activation in the pacing electrogram is filtered when this occurs. Fused beats occur when a paced beat and sinus beat occur close enough for their activations to overlap. Pacing subsequent to such sinus beats where they do not fuse are also discarded in the pacing electrogram, as the true activation occurs before intended and would therefore skew any conduction time value recorded for the associated pace. Entire pacing sequences are only discarded if the S1 stimulus preceding the extrastimuli, or the extrastimuli itself are not found; or if either are filtered in the pacing electrogram.
2.4. Isolate CT and EGM Window
Response activations in local electrograms are identified by isolating the largest peak in the CT window; where the CT window is defined by the interval between the preceding and subsequent stimulus from the pacing catheter (see Figure 6). The time delay between the stimulus (point A in Figure 6) and the largest peak in the resulting activation in the local electrogram (point B in Figure 6) defines CT. The CT window duration was restricted to an upper limit of 600 ms beyond the preceding stimulus. This was to account for instances when the subsequent stimulus is too far removed temporally, such as between the end of a pacing sequence and the start of the next one, and instances when a pacing beat did not capture.
An EGM window is defined as the duration between the center point between the preceding stimulus and response (point x in Figure 6), and the center point between response and subsequent stimulus (point y in Figure 6). This was done to minimize overlap between local activation and farfield activity due to pacing. Farfield activity refers to signal activity which occurs simultaneously to the stimulus, and does not represent local tissue activation; an example which can be seen on point C in Figure 6.
2.5. Feature Extraction
Within the previously defined EGM window, each peak in the local electrogram was assigned a score defined by the ratio of the integral of the peak divided by its distance from the center point; which denotes the largest peak in the window. The integral of the peaks was determined using the trapezoidal method  . Using this score, peaks were added to a range extending from the center of the window, until it encompassed 90% of the signal, as defined by the signal’s integral. This range denotes our ED. The assigning of scores to each peak was done using a divide and conquer recursion algorithm. CT is the time between stimulus and onset of local activation, which we define as the start of the ED. ND is determined by finding the number of extrema within the ED.
Due to their abundance, fractionation features obtained from the drive train were further refined using the following outlier removal equation using the modified Z-score, M, defined by Iglewicz and Hoaglin  :
Figure 6. Electrogram from roving catheter. Pacing spikes are shown by the clear triangles, while response spikes are shown by the red triangles. Sequence on the left side of the signal shows the steps involved in determining CT. The CT window is defined by the current and subsequent stimulus (A-C, shown by dotted line on left). Within this window the largest peak is found to isolate response spikes. Sequence on the right side of the signal shows the steps involved in determining the EGM duration. The EGM window (x-y, shown by the dotted line on the right) is defined by the half way point between the stimulus and response spike (D and E respectively, the location of which is identified when obtaining CT), and the half way point between the response spike and the subsequent stimulus (E and F respectively).
where defines the feature at index n, defines the median of sample, and MAD is the median absolute deviation. All values of M greater than 3.5 were removed as outliers, and the mean of the remainder of the remaining values was taken to represent the associated fractionation value for the train at the corresponding site. The difference between this value, and the CT determined for the extrastimuli was used to determine LT.
2.6. Physiological Constraints
Conventional mapping systems record the relative Cartesian coordinates of all points collected, which represent the geometry of the chamber being mapped. They also collect the relative coordinates of any catheter as it collects data. As such, a rough estimate can be made of the velocity between the pacing spike and the response spike by determining the distance on the shell between the RV and roving catheters during electrogram collection. To further improve accuracy, if the velocity was determined to be greater than 3 meters/sec (which corresponds to the fastest ventricular conduction physiologically observed  ), then the associated feature is discarded. Features are also discarded if, during data collection, the catheter is found to move drastically or if the catheter location is determined to be too far from the chamber geometry. If the algorithm is being processed in real time, as we would expect during an electrophysiological study, then any discards are flagged for user verification before either being removed (correctly identified discard), or reintroduced into the algorithm (incorrectly identified discard). If not in real time, then it is simply flagged for subsequent review. This also applies to any discards determined due the previously mentioned physiological interference.
2.7. Complexity Analysis
We assume that the length of the signal recorded is roughly consistent for all recorded sites and denote this length with n, and the total number of sites with m. The duration of the signal recorded will always be identical for all three electrogram channels used. Preprocessing for both pacingand roving electrograms can be summed as 6*O(n), as each Teagers energy operator and Butterworth filtering can be performed in O(n) operations. Isolating peaks in the pacing electrogram, and determining their corresponding response peaks in the roving electrogram, would involve one iteration of each signal. This would therefore have a complexity of 2*O(n). In determining activations in the surface ECG, a cross correlation would involve O(n*(n ? k + 1)) operations, where k is the length of our template; which, relative to the length of the signal is negligible and therefore can be reduced to O(n2). Finally, determining our ED would be an O(n*logn) process, as our recursive tree for determining scores splits into half upon each branch of each iteration (divide and conquer). Values associated with CT and ND are established in the extaraction of ED (onset of activation and number of extrema respectively); while determining LT is a simple arrhythmic process. Therefore, the sum of our complexity would be c*O(n2) +O(n*logn) +c = O(n2). Let c be a constant.
Data were collected from a total of 227 sites, 37.8 ± 8.8 locations of duration 32.1 ± 6.9 seconds per patient. A total of 264 electrograms collected from a roving catheter were manually and automatically annotated for fractionation. Of these, 60 were identified to have no redeemable features (physiologically induced noise, low contact, etc.) and we reseparated manually. All 60 were successfully discarded by our algorithm, yielding a specificity of 100%. Of the remaining 204 sequences, 16 were erroneously discarded by our algorithm, yielding a sensitivity of 92.16%. A comparison between annotations showed correlations of 0.98, 0.97, and 0.94 for CT, ED, and ND respectively. Subset comparisons for each pacing interval and feature, along with their respective r values (Pearson product-moment correlation coefficient), can be seen in Figure 7. A com- parison between the mean and standard deviations for these subsets can be seen in Table 1.
Figure 7. A comparison between the manually (y-axis) and automatically (x-axis) annotated data for 204pacing sequences.
Table 1. A comparison between all fractionation features over all cycle lengths; as annotated manually and automatically. Values are shown as mean ± standard deviation.
The meantime for five separate operators to manually annotate a pacing sequence for our fractionation features was 69.8 ± 19.8 seconds. The mean time it took for our algorithm to perform the same task over 1000 iterations was 0.15 ± 0.03 seconds. While it is to be expected that automation is faster, the point of interest to note is the magnitude of the discrepancy between the two values; which is substantial.
The use of automation in ventricular arrhythmias is not uncommon, and has been approached by several other investigators  -  . However, our method offers several key improvements. Many automation studies focus on detection of ventricular arrhythmias; primarily for their use in ICDs   . Thus, they are unable to detect arrhythmogenic sites, as ours does. Furthermore, by limiting the scope of our algorithm to a subset of substrate with ventricular arrhythmia, our implementation is optimized specifically to putative discrete areas of slow conduction resulting in re-entrant tachycardia.
While there have been many strides in the field, current cardiac models still cannot capture the stochastic nature of the electrophysiological conduction system of the heart  . Nor can animal models  . As such, validation through human data, as we have done here, yields the most correlative results for clinical situations. In humans, noise in the form of sinus and ectopic beats is not uncommon, and can lead to misleading premature activations. This can lead to incorrect conclusions about the electrophysiological makeup of the ventricle while mapping using this unique paradigm. As such, a single source of data for analysis is often inadequate. Furthermore, catheter movement during data collection can simulate fractionation, while low contact can mimic regions of low voltage  . Our use of both surface ECG, and cardiac and catheter geometry helps account for this and greatly reduces its interference.
The vast majority of automation studies focus on activation and substrate mapping   . But, as we previously discussed, recent studies showing efficacy of PEFA leads us to believe this might be a superior method for exposing and isolating fractionation a signature of slow conduction. Most techniques using the protocols similar to PEFA focus only on the feature extraction of CT as described originally by Saumarez et al.  . Our method expands on this by abbreviating the protocol and further extracting electrogram duration and deflections. Both of which have been independently shown to indicate possible ablation sites   . A collaboration of these features is likely to be much more indicative of arrhythmogenic sites than any single feature alone.
Observed results, for both manual and algorithmic annotation (Table 1), are consistent with previously observed trends associated with decrement pacing   . This leads credibility to both the use of our manual annotation as our gold standard, and to the validity of the output of our algorithm.
We have presented an automated fractionation detection algorithm which uses a modified PEFA protocol to identify the features: CT, ED and ND. Its performance has been shown to achieve a high degree of accuracy in comparison with manually annotated data, while greatly reducing associated annotation time. Furthermore, our automation revealed a high level of specificity and sensitivity. Current limitations in sensitivity may be overcome with operator intervenes, though future steps would involve refinement to further limit manual intervention. Nevertheless, we conclude that our approach provides a useful new method for automated detection of fractionated sites in patients with VT for use with commercial automated mapping systems, which can further be used to improve potential efficacy of ablation practices.