Fuzzy decision tree to classify complex fractionated atrial electrograms

Study population

The study cohort consisted of 11 patients who underwent catheter ablation of AF. The local Ethics Committee approved this study (according to the declaration of Helsinki), and all patients gave written informed consent. Eight patients had persistent and three paroxysmal AF. During the EP study, the left atrium was reconstructed using a 3D mapping system (Ensite NavX, St. Jude Medical, St. Paul, MN, USA) and PVI was performed. After PVI, bipolar electrograms were recorded using circular multipolar mapping catheters (10 polar Lasso, Biosense Webster, Diamond Bar, CA, USA; 14 polar OrbiterPV, Bard Electrophysiology, Lowell, MA, USA; 14 polar Optima, St. Jude Medical, St. Paul, MN, USA). After the procedure the electrograms and xyz-coordinates were exported and retrospectively analyzed. All together 605 recordings were analyzed. Each had a length of 5 s and was sampled at 1.2 kHz. The data was prefiltered by the measurement system with a high pass at 30 Hz and a low pass at 300 Hz.

For automatic classification, the Wells’ criteria [42] were modified by changing type 3 and 4. This led to a classification with a continuous increase in complexity (Figure 1).

Figure 1:

CFAE during atrial fibrillation.

Left: schematic presentation of the four variants according to the modified classification by Wells et al. Right: recorded electrograms during AF. C0 is a non-fractionated atrial electrogram of high frequency, C1 is a fractionated atrial electrogram with periodic activity, C2 is a mixture of periodic fractionated and periodic non-fractionated atrial electrograms, and C3 is a high frequency atrial electrogram with continuous activity.

Data preprocessing

All electrograms were preprocessed in the same way. Baseline wandering and low-frequency noise were extracted from the atrial electrograms using a discrete wavelet transform-based approach [16]. Based on the sample frequency (f_s) the signal is decomposed up to level n=ld(f_s) and the approximation of level n is set to zero. After this step, the signal is composed again. The high frequency disturbance is removed by a conventional Butterworth low pass filter of order 4 with a cut off frequency of f_LP=300 Hz [15]. After removing the noise by means of signal processing, the 605 stationary 5 s signals were classified by two physicians from different centers according to the above predefined classes (see Table 1). After the classification process, only electrograms with an unambiguous assignment were used to train and test the fuzzy decision tree (FDT). So finally the database consisted of 429 classified electrograms.

Descriptors of CFAE

To optimize the automatic electrogram characterization multiple features were generated. In the following sections a list of descriptors is delineated, reflecting the mathematical characteristics of CFAE.

Time domain descriptors based on non-linear energy operator: In sinus rhythm or other regularized atrial activity, intracardiac recordings are only presenting deflections if an excitation wave front propagates near the recording electrode. At baseline, no electrical activity is present. Therefore, the electrograms can be divided into segments with and without baseline crossings. Segments with baseline crossings were defined as active. Localization of active segments during AF is complicated, because they vary in length and number. A mathematical expression of active segments can be described as the energy, which forms the basis for time domain analysis.

Teager’s non-linear energy operator [NLEO, (1)] [38] forms the basis of this descriptor. Adding an adaptive threshold on the low pass filtered NLEO enabled the separation of active and inactive segments [28, 38]. The output of the NLEO E_n can be considered as an indication of the energy of the signal x(n) and it is proportional to the frequency and amplitude of the signal. Based on this segmentation a set of descriptors is defined.

(1)En=xn2-xn+1xn-1 (1)

The sum of the length of all active segments during a 5 s recording divided by the total length of the signal is called “activity ratio” (AR, descriptor D10).

(2)AR=1L∑i=1Mli (2)

with L the length of the total signal, M the number of active segments and l_i the length of an active segment i. Other descriptors are the mean length of active segments (MLAS, descriptor D11) and the standard deviation of active segments (sdMLAS, descriptor D12). They are defined as

(3)MLAS=1M∑i=1Mli sdMLAS=(1M∑i=1M(li-MLAS)2)1/2 (3)

Five more descriptors are defined based on the found active segments. These are the number of active segments (NoAS, descriptor D14), the mean number of zero crossings per active segment (ZCAS, descriptor D16), the standard deviation of zero crossings per active segment (sdZCAS, descriptor D18), the mean number of local maxima per active segment (LocMaxAS, descriptor D15) and the standard deviation of local maxima per active segment (sdMaxAS, descriptor D17). To analyze the signal curve within an active segment, a method is used that interprets the output of the NLEO as a probability density function. The standard deviation of this function is a measure of the concentration of the signal curve in time direction. The mean of the standard deviation of all active segments per signal is the descriptor mean variance in time direction (MVarTD, descriptor D1).

Phase space descriptors: A combined view of the signal’s magnitude and the alteration at a position t_i is enabled by the phase space. Here, the derivative of the signal x′(t_i) is plotted over its magnitude x(t_i). To weight the alteration and the magnitude equally, both are normalized to their maximum. To analyze the distribution of the samples in phase space, the phase space is divided into circular regions. Therefore a maximal distance d_max to the point of origin is defined. The distance is calculated by

(4)d(ti)=[x′(ti)2+x(ti)2]1/2 (4)

To be more robust against outliers d_max is defined as the mean of the 5% of the maximal distances (d_0.05max)

(5)dmax=10.05N∑j∈sd(tj), with  S={i|d(ti)∈d0.05 max} (5)

The aim of this definition is to use enough samples to lower the influence of outliers to the outer boundary, but with a minimum of samples. The phase space is divided into four regions depending on the distance d_max. Region 1 encloses all sample values within a distance of 0.05 d_max. The border for region 2 is at 0.1 d_max, for region 3 at 0.2 d_max, and for region 4 at d_max. The number of samples in each region divided by the total number of samples is called “Phase Space Sample Ratio” (PSSR 1, PSSR 2, PSSR 3, PSSR 4) (descriptors D3–D6). Therefore, the PSSR will group periodic electrograms (e.g., sinus rhythm) usually into region 1. For chaotic or non-periodic signals, the assignment is more random. To reflect this behavior a binary function is generated. This binary function equals “1”, if a sample is located in region 4 and “0” otherwise. As the entropy mirrors the information content of a random process, it is suitable to measure the periodicity of this binary function. The entropy of the binary function is termed as entropy of phase space (EPS 4, descriptor D7). To evaluate the periodicity of the binary function, the standard deviation of distances between two adjacent rising edges is calculated. The resulting descriptor is called mean cycle length of phase space (MCPS 4, descriptor D8).

Wavelet based descriptors: Fractionated activity exhibits complex and time varying morphology. This is reflected by the broad frequency spectrum present in those electrograms. Using a wavelet based approach allows the analysis of different frequencies on multiple scales.

According to the algorithm presented in [17] a descriptor is set up, that counts the zero crossings in a given signal as a measure for the fractionation of this signal. Therefore, the signal is decomposed into coefficients with wavelet Coiflet 4 up to level 10. As shown by [19] Coiflet 4 is a suitable wavelet to analyze fractionation of electrograms. The algorithm can detect regions with high magnitude and high slew rate. As intracardiac electrograms often differ in terms of frequency, the wavelet level with the largest magnitude is searched. According to the algorithm from [17] this level is multiplied with the two preceding (higher frequency) levels. As a result, active regions will be emphasized and inactive regions will be suppressed. In the resulting signal, x_res, fractionated regions will be detected by a search for zero crossings. The number of zero crossings in this wavelet based electrograms are used as a descriptor (FracSig, descriptor D13).

Similarity of active segments: Faes et al. [10] described an algorithm that compares the similarity of different regions of an intracardiac electrogram. Inspired by Faes’ work a similarity analysis is computed. The envelope of the absolute value of the analytical signal is calculated according to the algorithm described in [34]. The envelope follows the shape of the signal, but is always positive. This approach respects the increasing variation of the accurate shape in a more fractionated electrogram. Active segments are extracted using the NLEO-based segmentation algorithm. With these extracted segments a correlation matrix is built. First the segments are aligned using cross-correlation. On the overlapping parts, the absolute value of the correlation coefficient is calculated. On the basis of this correlation matrix, a clustering is performed. Starting with cluster of the two most similar segments, step by step the nearest similar segment is added and the similarity between the new segment and the cluster is computed. According to the definition of Kaufman and Rousseeuw [13] the similarity between two clusters is calculated from

(6)s (R, Q)=1|R||Q|∑i∈R, j∈Qs (i, j) (6)

where |R| and |Q| are the cardinalities of both clusters, s(i, j) is the similarity between element i from cluster R and element j from cluster Q. Here, the number of elements of cluster Q is set to 1. Finally, the mean value of the similarities of an active signal (Similarity AS, descriptor D9) is computed.

Amplitude statistics based descriptors: The histogram of an electrogram indicates the measured values during the time of recording. Little electrical activity in atrial electrograms is resulting in an amplitude around zero. So, the amplitude histogram has a high peak around “0”, whereas the amplitude distribution of electrograms with strong electrical activity is spread more widely and is more gaussian-like. This behavior can be described with the fourth standardized moment, the kurtosis [25]. To lower the impact of outliers the signal is divided into k segments with a length l_s. The global kurtosis (HistKurt, descriptor D2) is averaged over the kurtosis of the k segments. The segment length l_s is set to 1 s to make sure that even in sinus rhythm up to a heart frequency of 1 Hz there is at least one atrial excitation in this segment.

(7)HistKurt=1k∑j=0k-1kurtosis(x(t|j⋅ls≤t<(j+1)⋅ls)) (7)

Table 2 gives an overview of the 18 descriptors, most of them used by the FDT for classification of the electrograms.

Classification/fuzzy decision tree

Building the fuzzy decision tree: As a major advantage, FDT avoids sharp split nodes but instead assigns objects in a region around the split point a relative affiliation to the following child nodes. This fuzzy zone is realized using sigmoid functions of form:

(8)zl(ω)=1-11+exp{-σ(vω-s)} (8)

(9)zr(ω)=11+exp{-σ(vω-s)} (9)

as proposed by Chandra and Varghese [7]. z_l and z_r are the affiliation to the left and the right child node, s is the split point and v_ω is the descriptor value of descriptor ω. The coefficient σ is the standard deviation of this descriptor on the training data. If the descriptor values are scaled with a factor a, the width of the fuzzy zone is changing by a². This is not intended.

In this decision tree, we defined a coefficient k to replace σ. The beginning z_start of the fuzzy zone is set when the sigmoid function equals 0.01 and the end z_end is set when the sigmoid function equals 0.99. In general z_start (x)=g and z_end (x)=1–g. The range of a descriptor is the interval [x_min, x_max]. The width of the fuzzy zone is chosen as part p∈[0, 1] of the interval length i=x_max–x_min. It is assumed that the fuzzy zone with width p·i is set symmetrically to the split point s. This leads to

(10)zstart(s-p⋅i2)=g and zend(s+p⋅i2)=1-g. (10)

The sigmoid function is then expressed as

(11)z(x)=11+exp{-k(s-x)}. (11)

Inserting (10) in (11) and solving for k results in

(12)k1=-ln(g)-ln(1-g)-(p⋅i)/2, k2=ln(g)-ln(1-g)(p⋅i)/2, (12)

with k₁=k₂ one has

(13)k=-2[ln(1-g)-ln(g)p⋅i (13)

The boundaries for the interval [x_min, x_max] have to be determined from the training data. To be robust against outlier, the upper boundary value is set to 1.5 interquartile range of the upper quartile and the lower boundary to 1.5 interquartile range of the lower quartile. While training the decision tree, the aim is to find a split point at a node t, so that the relative membership of the objects to the right class is increasing. The Gini diversity index (GDI) is a regular criterion to evaluate possible split points [3]. The adaption of the GDI for usage in a FDT was made by Chandra et al. [7]:

(14)GDI(sj)=∑v=12N(tv)N(t)⋅[1-∑j=1J(Nj(tv)N(tv))2] (14)

where J is the number of classes; N(tv) is sum of the fuzzy-membership values of records of child node t_v with chosen split point s_j on descriptor j. N^(t) is the sum of the fuzzy-membership values of records in the tth partition before split. Nj(tv) is the sum of the product of fuzzy-membership values of the attribute and the fuzzy-membership values of the corresponding record for class j. The GDI is a rating criterion to find possible split candidates. To find those split points for a descriptor vector we do the following:

• sort descriptor values in descending order,

• set the positions where the class affiliation is changing as candidates for a split point,

• set the split point to the arithmetic mean of both neighboring descriptor values.

According to this algorithm for each descriptor vector the split point candidates are found. Finally, the candidate with the smallest GDI is chosen. The width of the fuzzy zone can be varied according to the distribution of classes over a descriptor. If there is a large overlap, the fuzzy zone should be larger than in the case of separated classes.

In this work, the best width of the fuzzy zone was determined by choosing the best result from a set of different fuzzy zones (Figure 2). For each fuzzy zone a 10×10 cross validation was realized and the mean error rate for each fuzzy zone was computed. The best error rate was achieved for a fuzzy zone with a width of 20%.

Figure 2:

Evaluation of different fuzzy zone widths to find the “best” fuzzy zone width. Displayed is the error rate as result of a 10×10 cross validation for each width.

To make a decision on the class belonging of a given record, the decision tree needs an inference instruction. In contrast to a decision tree with sharp split values, where a record will get the class membership of the leaf node to which it is assigned, here, a test record can be assigned to more than one leaf node by a FDT. The inference instruction for a given record X_Test considers the fuzzy membership to a leaf node z(t_i, X_Test) and its class membership c(t_i). The fuzzy membership of a leaf node results from the product of fuzzy memberships of the passed nodes. The total class membership of a record X_Test is defined as

(15)ctotal=(XTest)=∑ic(ti)⋅z(ti, XTest) (15)

The classification of a record results in partial affiliation to the CFAE classes. From these partial affiliations, the resulting CFAE class is derived by a majority decision. The classified record will get the CFAE class label of the class with the highest percentage. The percentage can also be shown as a measure of certainty.

Chandra and Varghese used different sharp stop criteria to limit the growth of the tree while training [7]. This technique has some limitations. Therefore, the stop criteria of this decision tree were given very weak definitions. This led to an oversized tree at first. In an additional step this oversized tree was optimized to the best size [3]. The aim of pruning was to optimize the tree in regard to the complexity and the predictive accuracy. A sequence of trees is produced by the use of cost-complexity-pruning [3]. Out of this sequence the best tree is chosen by evaluating the error rate in consideration with the size of the FDT. The final tree for the classification in a clinical context will be generated on the whole set of available classification data. The best tree is chosen again after cross validation and optimization by evaluating the size and error rate.

Validation and statistical analysis

The standard way to evaluate a classificator is to perform a 10×10 cross validation [3]. This means the data set has to be separated into 10 equal subsets. The tree is trained with nine subsets and tested with the remaining “unknown” subset. The correct rates and error rates are computed and rated. This “training, testing and evaluation” step has to be repeated until each of the 10 subsets was used to test the tree. This whole procedure is repeated 10 times with different randomly created subsets, which leads to a total of 100 trees that are evaluated.

To optimize this tree, the evaluation has not only been done on the CFAE database but also on datasets from the UCI Machine Learning Database [22]. The UCI database is a regular database to test and evaluate machine learning algorithms. From this database five data sets [Iris, Wine, Wisconsin Breast Cancer (WBC), Haberman’s Survival (HS), and Glas Identification (GI)] have been chosen to evaluate the FDT. Data sets with continuous feature values and multi-class assignments have been chosen to be most similar to the CFAE data. Finally the results have been compared to different decision tree algorithms (CART [3], SLIQ [23], C4.5 [35], FDTx [7]).

Evaluation of clinical data

To evaluate the distribution of the CFAE classes in clinical cases, bipolar electrograms were recorded for 5 s with a 300 Hz filter setting. Anatomical and electrical information were exported from the system and the 3D anatomical shells as well as the corresponding electrograms were reconstructed. At first, the algorithm was set to detect and display the mean cycle length (CL) of AF (mean dV/dt), with a deflection width of 10 ms and a refractory period of 37 ms. Electrograms with a mean CL>120 ms were color coded in violet, 120–70 ms in rainbow colors and <70 ms in white. The resulting map was compared to the CFAE classification map using the FDT.

Descriptors

All descriptors have been calculated using the database described previously. In the training phase, the FDT chooses the best descriptors to separate the classes depending on the cost function (e.g., GDI). The results of the 18 most frequently used descriptors are shown in Table 2.

For each CFAE class the median (μ˜) and the interquartile range (iqr) per descriptor is shown. The μ˜ and the iqr give an overview of the distribution of the descriptor values and the power to separate CFAE classes of each descriptor. Generally, there is no descriptor that can separate between all four CFAE classes alone. However, the descriptors can separate between different CFAE class groups. For example the descriptors PSSR 3 (D5) and AR (D10) can separate C0 from C3. The mean NoAS (D14) can be used to separate C0 and C3 from C1 and C2. The boxplots of these descriptors are depicted in Figure 3. To evaluate the results of the phase space descriptors, the numbers of samples per region are divided by the total number of samples per signal. The boxplot of the PSSR 3 (D5) is depicted in Figure 3. As can be seen, the median values are significantly different on a 5% level (boxplot notches) for all four CFAE classes. Also there is no overlapping of the iqr. The same applies to the AR. The range for the AR varies between 0 (no activity) and 1 (continuous activity). It increases from 0.18 (C0) to 0.32 (C1) and from 0.56 (C2) to 0.82 (C3). An increasing value is correlated therefore with an increase in entropy of the electrogram. The entropy of the number of hits in the different regions of the phase space (EPS 4) is an often used descriptor, too. EPS 4 can be used to distinguish between C0 and C2 and C3. The median values for NoAS (D14) are 18 (C0), 26 (C1), 24 (C2), and 15 (C3), respectively. The values are increasing from CFAE class C0 to class C1. The frequency and fractionation is also increasing from C0 to C1. For C1 and C2, the NoAS are nearly the same as these two states are similar with regard to the number of CFAE. For C3 the activity is getting more continuous and therefore the NoAS (D14) is decreasing again. In the case of continuous activity it will be one. Finally, the NoAS can separate C0 and C3 from C1 and C2. There is an overlap in the iqr of C1 and C2, but, although the value range is very similar the μ˜ are significantly different. For this four delineated descriptors the μ˜ are significantly different on a 5% level for all CFAE classes.

Figure 3:

Boxplots of four descriptors used by the decision tree depicted in Figure 4.

Red crosses mark outliers.

On the one hand, with increasing atrial activity the signals have more deflections and zero crossings and the iso-electric line is vanishing. Hence, the amplitude histogram is broad and the HistKurt (D2) has small or negative values. Sinus rhythm electrograms, on the other hand, mainly present an iso-electric line with some deflection other than zero. Therefore, the amplitude histogram has a pronounced peak. In this case the kurtosis values will be larger than for the ones for the continuous activity. This behavior is reflected by descriptors D1 and D2. When the heart rhythm evolves from sinus rhythm to AF the similarity of activation patterns gets lost [37]. As a consequence the active segments are getting dissimilar (D9).

Class C0 and C1 electrograms contain mostly regions without alteration (iso-electric line is prominent). This type of signals will be represented primarily in the phase space in region 1 and 2, whereas CFAE class C2 and C3 will be more prominent in region 3 and 4. Also the alterations of C2 and C3 signals make them lie more in the outer limits of the phase space. Descriptors D3–D8 mirror this. Going from the origin to the boundary of the phase space the number of C0 samples is decreasing while the number of C3 samples are increasing. PSSR 3 (region 3) contains 1% of C0, 4% of C1, 8% of C2, and 11% of C3 signals. EPS 4 and MCPS 4 give an information of the presence of a CFAE class in phase space region 4. These descriptors reflect the same outcome as the descriptor PSSR. The descriptors D13, D15, D16, D17, and D18, which reflect the fractionation of the signal, are increasing from C0 to C3, too. This corresponds to the increasing atrial activity and the increasing number of active fractionated segments (Figure 3).

Validation of the tree implementation

The FDT presented in this work was compared to the results of decision trees with sharp split values and the FDT by Chandra and Varghese (FDT_x) [7]. Table 3 gives an overview of the results. Presented is the error rate and standard deviation for each classifier after a 10×10 cross validation. The results for FDT_x and SLIQ were obtained from [7] and the results for CART and C4.5 were obtained from [31]. The decision trees CART, SLIQ, and C4.5 have sharp split values; the FDT by Chandra and Varghese uses a sigmoidal function as decision border. FDT is the tree presented in this work. In addition to the error rate, the fuzzy zone width for each data set is given. When comparing with other classifiers for the data sets Iris, Wine and WBC, the proposed FDT delivers comparable results. Focusing on the results of the Haberman’s Survival and the Glas Identification data sets, the FDT is significantly superior.

Resulting fuzzy decision tree

Applying a 10×10 cross validation on the presented FDT with a fuzzy zone width of 20% results in a mean correct rate of 80.65±3.3%. Correct rates for the different CFAE classes are C0: 83.1±4.4%, C1: 81.0±8.1%, C2: 75.8±8.4%, C3: 82.7±17.8%. The distribution of the wrong assignments of signals per class is shown in Table 4. Increasing the training data will lead to an increased accuracy of the FDT. The depicted tree in Figure 4 is chosen from the sequence of cross validated trees. This tree has a correct rate of 86.1% (Figure 4).

Figure 4:

Example of the optimal tree chosen from the cross validation process.

At each node the chosen descriptor and its split value is depicted. For each leaf node the membership result is shown. The correct rate for this tree is 86.1%.

Visualization of CFAE classes

Using the CFAE mean algorithm, areas with high atrial frequencies were located in different areas of the left atrium. A loose correlation between CFAE class 2 and 3 could be estimated. Interestingly, in the FDT map electrograms with continuous activity (C3) were surrounded by C2 and C1. C3 were present mainly around the PV ostia, the anterior wall and the inferior part of the posterior wall (Figure 5).

Figure 5:

Example of a clinical evaluation of the CFAE classes using the FDT in comparison to a CFAE mean map.

View of the posterior wall with left upper and lower pulmonary vein (left) and right upper and lower pulmonary vein (PV) (right). (A) automated CFAE mapping algorithm (Ensite NavX). Bipolar electrograms were recorded for 5 s, filter settings were at 30–500 Hz. The algorithm was set to detect and display the mean CL of atrial fibrillation (mean dV/dt), with a deflection width of 10 ms and a refractory period of 37 ms. Electrograms with a mean CL>120 ms were color coded in violet, 120–70 ms in rainbow colors and <70 ms in white. (B) CFAE classes using the FDT. C0 (blue): non-fractionated atrial electrogram of high frequency, C1 (green): fractionated atrial electrogram with periodic activity, C2 (yellow): mixture of periodic fractionated and periodic non-fractionated atrial electrograms; C3 (red): electrogram with continuous activity. Electrograms with continuous activity C3 are surrounded by C2 and C1. C3 were present mainly around the PV ostia and the inferior part of the posterior wall.