Complexity of Atrial Fibrillation Electrograms Through Nonlinear Signal Analysis: In Silico Approach

Abstract

Identification of atrial fibrillation (AF) mechanisms could improve the rate of ablation success. However, the incomplete understanding of those mechanisms makes difficult the decision of targeting sites for ablation. This work is focused on the importance of EGM analysis for detecting and modulating rotors to guide ablation procedures and improve its outcomes. Virtual atrial models are used to show how nonlinear measures can be used to generate electroanatomical maps to detect critical sites in AF. A description of the atrial cell mathematical models, and the procedure of coupling them within two‐dimensional and three‐dimensional virtual atrial models in order to simulate arrhythmogenic mechanisms, is given. Mathematical modeling of unipolar and bipolar electrogramas (EGM) is introduced. It follows a discussion of EGM signal processing. Nonlinear descriptors, such as approximate entropy and multifractal analysis, are used to study the dynamical behavior of EGM signals, which are not well described by a linear law. Our results evince that nonlinear analysis of EGM can provide information about the dynamics of rotors and other mechanisms of AF. Furthermore, these fibrillatory patterns can be simulated using virtual models. The combination of features using machine learning tools can be used for identifying arrhythmogenic sources of AF.

1. Introduction

The most common sustained cardiac arrhythmias in humans are associated with the atria. Atrial arrhythmias, mainly atrial fibrillation (AF), frequently provoke incapacitating symptoms and severe complications such as stroke and heart failure [1]. Overall, 20–25% of all strokes are caused by AF [2]. The presence of AF is related to a significant increase in morbidity and mortality [3]. Electrocardiogram‐based surveys suggest that 1% of the total population is affected by AF [4].

There are a large number of clinical conditions that are associated with an increased incidence of AF. This contributes to a progressive process of atrial remodeling characterized by a set of changes in atrial properties that contributes in sustaining of the arrhythmia. These changes include alterations in the electrical cellular activity, calcium handling and in the atrial structure such as cellular hypertrophy and fibrosis. They have been described in some animal models [5–8] and in humans [9–11]. These alterations may favor the occurrence of triggers that initiate the AF and the formation of a substrate that promotes its perpetuation. Changes in electrical activity cause a significant shortening of the action potential duration (APD) and a decrease in refractoriness [8–10], which may support the initiation and maintenance of multiple re‐entrant waves, as suggested by experimental studies [5, 9].

It is well known that AF can be caused by different mechanisms, including single‐circuit re‐entry, multiple‐circuit re‐entry, rapid local ectopic activity and rotors [12–15]. It is very important to know the mechanisms underlying AF, since these have implications in the treatment of the disease. An important percentage of patients suffers of paroxysmal AF, which is initiated by focal triggers that are localized at preferential sites, mainly in the pulmonary veins (PV) [13]. Electrical isolation of pulmonary veins can prevent recurrence of AF in 70–80% of these lone AF patients. The rationale for this is the crucial observation, reported in [13], that AF was mostly triggered by ectopic beats arising from the muscle sleeves of the pulmonary veins. They demonstrated that atrial rapid paces or ectopic activity originated in the proximities or in the interior of the pulmonary veins could act like triggers, and, in some cases, they would be responsible for the maintenance of paroxysmal AF episodes [16, 17]. A unifying theory suggests that rapid focal activity is responsible for generating atrial, which is necessary to maintain a substrate for the generation of multiple re‐entrant waves [18, 19]. While paroxysmal AF is maintained predominantly by ectopic focal activity or local re‐entrant circuits located in one or more pulmonary veins, as the arrhythmia evolves into more persistent forms promoted by atrial remodeling, the mechanisms that maintain AF move toward the atria and are increasingly based on re‐entry substrates [11, 20–22]. Based on clinical [23–25] and experimental [14, 26] results, certain types of AF can be attributed to a stable high‐frequency rotor or a small number of rotor waves in left atrium, which maintain the arrhythmia, whose periodic activation can be converted into a chaotic pattern when the wavefronts propagate across the atrial wall. This phenomenon, known as the mother rotor hypothesis, is the most recently proposed mechanism of AF [27], which suggests that AF is triggered by a series of ectopic beats, whose wave fronts give rise to a rotor. The rotor is a stable re‐entry around a functionally unexcitable core [15] that works as a maintenance mechanism with some spatial temporal stability, activating the local tissue at high frequency, generating wave fronts that fragment and propagate as multiple daughter wavelets. Stable rotors are at diverse locations, mostly in the left atrium, including sites outside the pulmonary veins, as well as the posterior, inferior, and roof regions. Several studies have observed rotors in in vitro and animal models [14, 28, 29], and its presence in humans has been reported [27, 30, 31].

All localized sources of AF as re‐entry circuits, ectopic foci, or rotors, cause fibrillatory conduction remote from the source, which is difficult to distinguish from propagation that maintains the AF by multiple wavelets, and any of these phenomena may generate rotors registered by intracardiac recordings [32, 24]. Studies carried out by [33, 34] have shown evidence that human AF can be sustained by localized rotors; however, the existence of stable rotors maintaining AF remains a focus of discussion. Hence, the importance of the implementation of virtual models and computational tools in the studies of AF sources and its relation with the signals recorded from surface. In the same way, the rotors have been shown in computer models of AF [35–38].

Ablation has emerged as an important treatment strategy for AF [39]. Pulmonary veins isolation reaches single procedure success rates of 60–80% for paroxysmal AF treatment [40]. However, for the chronic case, this strategy does not achieve satisfactory outcomes [41]; consequently, complex ablation lines are added to the procedure [45]. Catheter ablation guided by electroanatomic mapping has revolutionized the treatment of permanent AF. Identification of AF mechanisms could improve the rate of ablation success. However, the incomplete understanding of those mechanisms makes difficult the decision of targeting sites for ablation. To overcome this limitation, two‐dimensional and realistic three‐dimensional human atria computer models have been developed to investigate the relationship between the characteristics of EGM and the propagation pattern associated with them [42–45].

It is thought that the different mechanisms lead to changes in the characteristics of spatiotemporal organization of AF [44, 45]. EGM‐guided ablation procedures have been proposed as an alternative strategy, which involves mapping and ablating focal sources or complex fractionated atrial electrograms (CFAE) [46]. Recent studies have shown that a large number of sites representing AF substrates are characterized by a high degree of disorganization in EGM signals [29], and therefore, signal processing methods are being designed in order to quantify their degree of fractionation [47, 48]. The relationship between the rotor tip and CFAE has been published in recent studies [36, 45, 47, 49, 50]; however, automatic rotor mapping methods have not been fully developed.

Different mapping techniques are being used to identify target sites for ablation, and some of them are activation waves, voltage, dominant frequency [23, 51–55] and CFAE maps [39, 46, 56]. Nevertheless, there are several limitations with these techniques, and one of the most important is that they depend considerably on the electrophysiologist expertise [57].

The term CFAE was introduced by Nademanee [46] as a pathophysiological concept; however, its definition is unclear and broad and includes inherent subjectivity [58]. CFAE are formally defined as follows: (1) atrial EGM that have fractionated electrograms composed of two deflections or more, and/or perturbation of the baseline with continuous deflection of a prolonged activation complex over a 10 s recording period; (2) atrial EGM with a very short cycle length (<120 ms) over a 10 s recording period. CFAE definition does not distinguish between different morphologies; consequently, fractionated EGM according with CFAE definition are not necessarily related to arrhythmogenic substrates. Moreover, several studies have shown conflicting results with Nademanee [56, 59]. This may lead the electrophysiologist to confuse the fractionated EGM that are functional in nature [60] with fractionated EGM corresponding to arrhythmogenic sources, leading to incorrect target sites for ablation. This also makes clinical results difficult to compare. Inconsistent results have been found in different studies using the CFAE concept. Although the concept of CFAE has made a great contribution to the study of AF, it may fail to describe the wide range of fractionation that occurs in the EGM signals in specific cases. Thus, the electrogram‐guided approach by mapping and targeting areas of CFAE is a technique that is still debated [61], and further studies are needed to strengthen the treatment techniques.

To date, indexes computation from EGM is carried out mostly based on the detection of local activation waves and time intervals. However, the concept of EGM fractionation, such as CFAE, defined using time criteria, is not enough to describe the complexity of EGM during AF. Morphological irregularity and temporal variability of the signal must be included. To overcome the limitations of CFAE definition, designation of different degree of fractionation has been proposed based on the perturbation of baseline and the presence of continuous deflection [36]. Some studies have presented evidence linking the vortex of the rotor with a high degree of fractionation in the EGM [36, 47]. However, the patterns of EGM signals associated with the rotors are still unknown, and the development of tools for characterizing fractionated EGM signals is still a current topic. In addition, there is not a standard for description of fractionation degrees, and there are different descriptions in the literature. In order to quantify the behavior of different morphologies of EGM signals and the fractionation degrees, recent studies have helped to understand the concept of EGM fractionation as a nonlinear phenomenon [47, 48, 62]. Mathematical descriptors of the nonlinear dynamical behavior have been used to study cardiac signals and are based on the description of the state of the system and its evolution. If a linear law does not describe the evolution of the system, the dynamics is nonlinear. Thus, indexes and features calculated using nonlinear dynamics theory become an alternative for enhancing assessment of EGM complexity.

2. Simulated atrial fibrillation dynamics

The mathematical modeling of atrial electrophysiology has become a useful method for studying the underlying mechanisms responsible for AF. Different models of human atrial electrophysiology have been published with various formulations of ionic currents and calcium handling, therefore, with different electrophysiological properties [63–67]. In our studies, the transmembrane potential was based on the Courtemanche‐Ramirez‐Nattel and Kneller [63, 68] model of human atrial cell kinetics in the presence of 0.005 μM of acetylcholine (ACh). Based on the experimental data [9, 10, 69], the cell model was modified in order to reproduce electrophysiological conditions of permanent AF: the maximum conductance of delayed rectifier potassium current (I_Kur) and transient potassium current (I_to) was decreased by 50%, the maximum conductance of potassium time independent current (I_K1) was increased by 100%, and the maximum conductance of L‐type calcium current (I_CaL) was decreased by 70% (see Table 1). Thus, the action potential duration at 90% of the repolarization (APD90) was reduced by 70% (Figure 1A), which is in accordance with experimental studies developed by Workman et al. [9] and Bosch et al. [10] in isolated myocytes from patients with permanent AF, using the whole cell patch clamp technique.

Figure 1.

(A) Action potential under physiological and permanent AF conditions. (B) S1–S2 cross‐field protocol. The plane S1 and rectangular S2 stimulus are shown. (C) Stable rotor under permanent AF conditions. (D) Contour map, the circular point indicates the rotor tip location.

Conductance	Permanent AF
gIKL	Increased by 100%
gIKACh	–
gIKuv	Decreased by 50%
gIto	Decreased by 50%
gICaL	Decreased by 70%

Table 1.

Conductance (g) changes of currents for permanent AF.

The development of multi‐dimensional models has allowed representing fibers, pieces of tissue and atrial cavities. These models have served to approach, among other aspects, the biophysics of the action potential propagation under physiological and pathological conditions.

Different studies have developed and implemented two‐dimensional (2D) models of atrial tissue to study the action potential propagation dynamics under physiological and pathological conditions [36, 37, 68, 70].

We have developed a two‐dimensional model of human atrial tissue. The tissue surface was discretized into a 150 × 150 hexahedral mesh (22,500 elements and 45,602 nodes). Spatial resolution was 0.4 mm. The electrical propagation of the atrial action potential was modeled using the monodomain reaction–diffusion equation:

$\frac{\text{1}}{S_v}\nabla \cdot (D\nabla V_m)=C_m\frac{\mathrm{\partial }{V}_m}{\mathrm{\partial }t}+I_{ion}-I_{stim}$E1

where S_v corresponds to the surface‐to‐volume ratio, D is the conductivity tensor, C_m is the specific membrane capacitance (50 pF), I_ion is the total ionic current that crosses the membrane cells, V_m is the membrane potential, and I_stim is the stimulus current. The electrophysiological model was integrated into the two‐dimensional virtual model. The equation was solved using the EMOS software, which is a parallel code that implements the finite element method and operator splitting [42, 71] for solving the monodomain model of electrical propagation and allows to calculate EGM as postprocess. The time step was fixed to 0.001 ms.

The tissue was considered isotropic. A conductivity of 0.3 S/cm was assigned to obtain a realistic conduction velocity of 60 cm/s.

Rotor was generated by S1–S2 cross‐field stimulation protocol. Stimuli pulses were rectangular with 2 ms of duration and 6 mA of amplitude. The S1 was a train of five plane stimuli applied at the left boundary of the model at a basic cycle length of 1000 ms. The S2 stimulus was rectangular (3 × 2 cm) and was applied 40 ms after the last S1 in the inferior left corner of the model (Figure 1B).

A 4s stable clockwise spiral wave (rotor) was observed in the two‐dimensional model after applying the protocol (Figure 1C). The APD shortening due to the atrial remodeling allowed the rotor stability over time. Wijffels et al. [8] in goats demonstrated that maintenance of AF by pacing in the normal goat heart resulted in the development of sustained AF within 1–3 weeks. This observation of tachycardia‐induced electrical remodeling creating a background for persistent AF led to the concept that “Atrial Fibrillation Begets Atrial Fibrillation.” The longer duration and stability of the AF episodes was explained by a shortening of the wavelength of the atrial impulse. The simulated rotor pattern agrees with the “rotor hypothesis” proposed by Jalife [15]. The contour map was implemented to show the core of the rotor (circular point in Figure 1D).

Virtual models of cardiac structures are needed in the study of atrial arrhythmias that are critically dependent on the spatial organization of cardiac structures and fibers. Some atrial arrhythmias, such as AF, need complete structures to perpetuate them. The validity and enormous potential of multi‐scale heart models in improving diagnosis, prevention and therapy of cardiac pathologies are supported by different scientific results. These virtual models have opened new horizons in the study of the complex mechanisms underlying atrial arrhythmias and their treatment, either pharmacological or surgical. The earliest complex atrial models incorporated parts of the anatomical architecture of the atria [72, 73]. During the last few years, a significant number of new models of the animal and human atrial anatomy have been published [44, 74–88]. Some of these models included the atrial fiber architecture observed in histological sections. Studies developed by Seemann et al. [75] and Aslanidi et al. [79] included only the main bundles fiber orientation in their human atrial model. An image‐based anatomical model of the sheep atrial has been published [80] including realistic fiber orientation. The model reproduces the whole atria with highly detailed myofiber architecture. The three‐dimensional models also allow relating the arrhythmic behaviors as focal activity, rotors and multiple wavelet reentries, with their manifestation in the electrograms [44, 81].

A realistic three‐dimensional model of human atria including the main anatomical structures (Figure 2A), electrophysiological heterogeneity, anisotropy (conduction velocity in the direction of myocardial fibers is usually several times larger than that vertical to them), and fiber orientation (Figure 2B) was developed in an earlier work [42]. It includes 52,906 hexahedral elements (polyhedrons with six faces, eight corners or nodes, topologically equivalent to cubes). The mathematical atrial cell model coupled within three‐dimensional (3D) virtual atria model was used to simulate AF dynamics. AF episodes were generated by the S1–S2 stimulation protocol as follows [36]: a train of five stimuli with a basic cycle length of 1000 ms was applied in the sinus node area for a period of 5 s to simulate the atrial sinus rhythm (S1). Based on the study developed by Haissaguerre et al. [13], a burst pacing of 6 ectopic beats to high frequency (S2) at cycle length (CL) of 130 ms was delivered into the right superior pulmonary vein after the last S1.

Figure 2.

(A) Frontal view of the three‐dimensional model of human atria. (B) Fiber orientation. (C) AP for different atrial areas (CT: crista terminalis, PM: pectinate muscles, APG: appendages, AVR: atrioventricular rings, and AWM: atrial working myocardium) under physiological conditions. (D) Activation isochronal maps. Stable rotors located in the posterior wall of the left atrium (#2) and A in the superior vena cava (#1) are showed. A block line can be seen at the right inferior pulmonary vein (#3).

During AF activity initiated by the ectopic activity, it was observed the generation of two rotors of stable activity during 5 s of simulation. One was located in the posterior wall of the left atrium, near the left pulmonary vein (#2 in Figure 2D), and the other was located in the superior vena cava (#1 in Figure 2D). The rotors were generated spontaneously, which is in agreement with the rotor hypothesis [27]. Additionally, a block line located over the inferior right pulmonary vein has been observed (#3 in Figure 2D).

3. Simulated atrial electrograms

Several studies [21, 60, 82–84] have investigated the effects of factors such as slow conduction, anisotropy, conduction blocks, re‐entries, and wave collisions, on the morphology of unipolar and bipolar EGM. However, it is still not entirely clear to what extent these factors contribute to temporal and spatial variations in EGM morphology as observed during AF.

Calculating atrial EGM from the virtual models allows the study of EGM morphology and their relationship with arrhythmogenic sources.

Unipolar EGM are modeled as the register of the extracellular potential measured by a positive polarity electrode whose reference (zero potential) is located at infinity. The distance from the electrode to the surface quantifies the influence area of the electrode, so the closer it is to the tissue, the greater the field uptake. The extracellular potential (Ф_e) was computed using the large volume conductor approximation [85, 86]:

${\phi }_e\left(r\right)=-K\mathrm{ʃ}\mathrm{ʃ}\mathrm{ʃ}{\nabla }^{\prime }{V}_m\left(r^{\prime }\right)\cdot {\nabla }^{\prime }\left[\frac{1}{r^{\prime }-r}\right]dv$E2

where K is a constant that includes the ratio of intracellular and extracellular conductivities (${\sigma }_i/4\pi {\sigma }_e$), ∇′V_m is the spatial gradient of transmembrane potential V_m, r is the distance from the source point (x, y, z) to the measuring point (x′, y′, z′) and dv is the differential volume.

Zlochiver et al. [45] investigated the regularity of EGM in the presence of stable rotors, in a two‐dimensional atrial model. Jacquemet et al. [87] in a computer model representing a monolayer of atrial cells concluded that microscale obstacles cause significant changes to EGM waveforms. Using two‐dimensional computer models and cell cultures, Navoret et al. [88] detected CFAE using the criteria of cycle length, number of deflections, and amplitude. They established a relationship between the detected CFAE and the presence of rotors and shock waves, but they failed to differentiate them. Ashihara et al. [89], using a two‐dimensional myocardial sheet of size 4.5 × 4.5 cm, studied the role of fibroblasts in CFAE during AF. Yun et al. [90] reported that CFAE in a homogeneous two‐dimensional AF model were weakly correlated with wave break, phase singularity, and local dominant frequency.

We simulated that in the two‐dimensional atrial model, a total of 22,500 virtual electrodes (150 × 150, one for each element of the model), spaced by 0.4 mm at a distance of 0.2 mm above the atrial surface, unipolar EGM were calculated with temporal resolution of 1 ms.

The 98.9% of EGM, located away from the rotor tip, present simple morphology (Figure 3A). The remaining 1.1% of EGM, located at the rotor tip, exhibits potentials composed by two or more deflections (Figure 3B).

Figure 3.

(A) Regular EGM calculated in ‘*’ from Figure 1C. (B) Fractionated EGM calculated in the rotor tip.

The mechanism by which fractionation of unipolar EGM occurs in our simulations can be explained as follows: the rotor is a singularity point or phase singularity, when the rotor is stable it pivots around a circular trajectory forming the core of the spiral wave, afterwards the pivot point is affected by the wavefronts from the rotor tip. When the wavefront passes near to the pivot point in each rotation cycle, several electrotonic potentials (nonpropagated local potential) are observed; consequently, irregularity and fractionation arise [36]. Our results are consistent with other studies, in which unipolar EGM symmetry was affected by the wavefront curvature (convex, concave or amorphous) [44] and fractionated unipolar EGM were observed at pivot points (functionally unexcitable core around which the rotor turns) [82]. Umapathy et al. [50] reported that CFAE were located in the region of a rotor tip and sites where wave breaks, using a murine HL‐1 atrial monolayer model.

Most of the in silico studies using three‐dimensional atrial models have characterized the simulated arrhythmias by observing the re‐entrant patterns. Few authors [44] have also calculated EGM in a circular region on the free wall of the right atrium, using the 16 unipolar virtual electrodes on a simplified three‐dimensional model of human atria. They suggested that analysis of the amplitude and symmetry of unipolar atrial electrograms can provide information about the electrophysiological substrate maintaining AF. Hwang et al. [81] calculated bipolar EGM in a personalized three‐dimensional left atrial model in order to applied virtual ablation at CFAE points.

We calculated 42,835 EGM in the whole atrial surface of the three‐dimensional atrial model, over a 4‐s window and recorded at 1 kHz. Bipolar EGM were calculated by subtracting two 1 mm‐spaced adjacent unipolar EGM.

Fractionated atrial EGM were shown to be located in rotor tip areas, when the tip of the rotor turned on this point, displaying low voltage and irregular morphology with potentials composed by two or more deflections (Figure 4A and B). The wavefront of the rotor surrounds the pivot point, without depolarizing it completely, which results in multiple low amplitude deflections in the EGM.

Figure 4.

EGM calculated in the three‐dimensional atrial model, under simulated AF episode. Fractionated EGM corresponding to the rotor in the posterior wall of the left atrium (A) and in the vena cava (B). (C) EGM corresponding to a functional block line. (D) Regular EGM corresponding to plain wavefront.

The EGM corresponding to the block line present fractionation; however, the activation patterns are visible, and their amplitudes are similar to nonfractionated EGM (Figure 4C). The EGM from sites with a plain wavefront are regular with potentials composed by one deflection (Figure 4D).

We identified the area in the posterior wall of the left atrium where the rotor spins (shaded circle in Figure 5A). EGM signals obtained from this area were used. From the selected region of the model, a conversion was made from the three‐dimensional coordinate system to the two‐dimensional coordinate system (x, z), taking advantage of the very low dispersion in y. EGM were converted to bipolar EGM, and this task was accomplished by creating a virtual mesh with 1 mm spacing, and performing a match with the two‐dimensional model surface. In this way, the difference between two adjacent signals in the mesh was calculated, obtaining a bipolar signal. In the same way, the results show that the rotor vortex area is associated with signals presenting high degree of fractionation (Figure 5B). For the contrary, the EGM from sites with a plain wavefront are regular with potentials composed by one deflection, similarly to the unipolar EGM morphology (Figure 5C).

Figure 5.

(A) The area in the posterior wall in the model to obtain bipolar EGM. Samples of two bipolar EGM are shown, a fractionated signal from the rotor tip (B) and a regular activation pattern from nonrotor area (C).

4. Estimation of nonlinear features for electroanatomical mapping

EGM‐guided ablation has been proposed as a strategy to find critical sites of AF as target sites for ablation. Multiple clinical trials have shown that ablation of fractionated electrograms adds no benefit to conventional AF ablation with pulmonary vein isolation [91, 92]. This is likely because sites of electrogram fractionation, according with CFAE definition, not always correlate with sites of arrhythmic drivers and can also represent sites of wavefront collision or slow conduction, among others. Although CFAE may be relevant to detect areas that maintain AF, further characteristics apart from fractionation should be important to identify the atrial sites that maintain AF.

To overcome the limitation of CFAE, nonlinear analysis of EGM signals has been proposed by several authors to analyze the signals using further characteristics apart from time intervals or number of deflection [93, 94]. Nonlinear features are studied using the raw EGM signals, and it is not necessary to detect local activation waves. This is an important property, because in fragmented signals detection of activation waves is not always feasible. Nonlinear features as entropy estimation and fractal analysis are compute over each single EGM, and its value is related to the complexity of the signal.

Nonlinear mathematical tools can be used to quantify the irregularity of a signal. During the last years, measures have been developed to estimate the complexity of biomedical signals. Main goal of these advancements is to provide new theories about the dynamics of biological systems. Such is the case of Kolmogorov‐Sinai entropy, Lempel‐Ziv complexity, and correlation dimension, among others [95]. However, most of the complexity indexes require long time series to obtain reliable and convergent measures. Pincus [96] proposed the statistic approximate entropy (ApEn), which solves the issue of short time series, and it is aimed to measure the complexity degree and the presence of similarity patterns. Further developments have followed the ApEn, taking this as a starting point, such as the sample entropy [97], the fuzzy entropy [98] or hierarchical entropy [99]. Some authors have reported the use of nonlinear features to evaluate their suitable for locating critical sites in AF. For instances, Ganesas et al. [47] reported that sites near to the rotor tip present high values of Shannon Entropy (ShEn) in EGM signals recorded from cell cultures and simulated episodes of fibrillatory conduction in two‐dimensional models of atrial tissue.

4.1. Approximate and Shannon entropy definitions

In general words, entropy has been conceived as a measure of the degree of disorganization or irregularity of a process. The most organized the process is, the lower the entropy related to it.

The statistic $ApEn(m,r,N)$ depends on the length N of the time series $x\left(n\right)$ (where n is), the positive integer m (where $m\le N$) and the positive real number r. Defining:

${\Phi }^m\left(r\right)=\frac{{\sum }_{i=\text{1}}^{N-m+\text{1}}\text{log[}{C}_i^m\left(r\right)\text{]}}{N-(m-\text{1})}$E3

we have that $ApEn(m,r,N)={\Phi }^m\left(r\right)-{\Phi }^{m+\text{1}}\left(r\right).$

The variable $C_i^m\left(r\right)$ counts the number of segments of length m that are within the boundaries defined by r. Thus, ApEn(m, r, N) measures the logarithmic frequency of the tool measures the logarithmic frequency that those segments of length m that are close remain close after increasing the length of the segments by one. In such a way, the statistic ApEn provides a measure of irregularity of the signal, implying strong regularity when ApEn value is small, and irregularity when ApEn value is large [96].

In previous work, we have reported the use of ApEn to evaluate the location of rotors and block lines in a three‐dimensional model of human atrial [100]; and the use of multifractal analysis as a tool to discriminate between four levels of fractionation according with a modified Well’s approach [101]. In this work, we tested several nonlinear features using EGM signals recorded from a two‐dimensional model of atrial tissue and a three‐dimensional model of human atrial. Additional, we test the used of combination of nonlinear features using clustering method to study the distribution of different EGM patterns over the atrial surface.

Another index that estimates the entropy value from an N‐point signal $x\left(n\right)$ is Shannon entropy (ShE) defined as:

$ShEn=-\sum _{i=1}^N{p}_i{\log }_2\left(p_i\right)$

where p_i is the probability of assuming the corresponding x(i) value. Both, ApEn and ShEn, consider that a high value of repeated patterns implies order. Thus, they make their respective estimations of a signal irregularity by counting repetitive patterns, where the ApEn has a more elaborate method of defining and counting these patterns.

4.2. Nonlinear features to complexity estimation for building maps: two‐dimensional model case

Nonlinear features such as ShEn and ApEn have been used and tested for locating stables rotors in simulated episodes of fibrillatory conduction. We have tested ShEn maps and ApEn maps using the EGM signals from the two‐dimensional model. Results of this approach are shown in Figure 6A and E. These maps are constructed with high resolution using all the signals available in the model: 22,500 with spatial resolution of 0.4 mm. However, in the real case, the resolution of the electrodes can be lower, so in Ref. [102] was carried out a study about the EGM maps analysis reducing their resolution. Figure 6 shows the maps of two‐dimensional model of AF reconstructed from the entire model with a 75% of reduction of the electrodes number (resolution: 37 × 38) and characterized using the features ShEn and ApEn of the EGM signals, respectively. A reconstruction of the entire model was developed using the interpolation techniques: Inverse distance weighted ‐IDW [103] (Figure 6B and F), IDW with Mean Filter–MF [104] (Figure 6C and G), and backpropagation artificial neural network—BPANN (Figure 6D and H). The best result is obtained with BPANN algorithm. Backpropagation artificial neural network (BPANN) is a type of artificial neural network that assumes the function of a common and complex nervous system, and BPANN is widely used in machine learning for clinical research [105, 106]. BPANN is trained using Levenberg‐Marquardt backpropagation algorithm [107]. This technique was applied for predicting the values of unknown points in order to increase the resolution of the two‐dimensional Map. BPANN has a structure (layers and neurons—[2 5 4 3 2 1]), which was defined applying a heuristic adjustment based on the minimum error for the mean map. The performance was assessed using the root mean squared error (RMSE).

Figure 6.

(A) Two‐dimensional ShEn map. ShEn map reconstructed 75% using IDW‐MF (B) and BPANN (C). (D) Log entropy map. (E) Two-dimensional ApEn map. ApEn map reconstructed 75% using IDW-MF (F) and BPANN (G). (H) Log entropy map.

5. Approximate entropy for rotor detection: three‐dimensional model case

Motivated by the features of the ApEn, as a signal analyzing tool, and its important presence in several studies of complex biological systems [108–113]; our group has performed a study that relates AF mechanisms, such as rotors, with high degree of irregularity of EGM by means of the ApEn. In the following sections, ApEn theoretical definition and its interpretation are presented. Moreover, ApEn electroanatomical maps obtained from the virtual models are presented, as well as the feasibility of characterize fibrillatory mechanisms in space and time. It follows a detailed analysis of our results and their implications.

5.1. Approximate entropy for locating critical sites in AF

As stated earlier, our group has performed numerical experiments to assess the regularity of atrial EGM signals by means of the ApEn. Our research is based on three hypotheses: (1) Fractionation of EGM increases the ApEn values. (2) High ApEn values can be related to the tip of a rotor. (3) Information about spatial and temporal dynamics of a rotor could be obtained using moving window ApEn.

In order to calculate the ApEn values from the virtual unipolar EGM, we define the parameters m = 2 and r = 0.1 according to the interval of values suggested by Pincus [96], and N = 1000. Figure 1 shows three EGM of 1000 points each, corresponding to minimum, intermediate and maximum ApEn values. The ApEn corresponds with the morphological complexity of the EGM: High values of ApEn mean irregularity or fragmentation of the EGM, and vice versa. In Figure 7, the EGM of the bottom present fragmentation of activation waves and baseline irregularity. Intermediate values of ApEn represent transitions between nonfragmented and fragmented EGM, in which differences between the patterns of activation waves can be seen, as can be seen in the EGM of the middle.

Figure 7.

Three degrees of EGM irregularity and the corresponding ApEn value.

Figure 8 (middle) shows the electroanatomical map of ApEn(2, 0.1, 1000) for the first second of the episode. The areas of R1 and R2, right inferior pulmonary vein and coronary sinus, have high ApEn values (red). The ApEn values (green) increase in the appendix, in pulmonary veins, and on the posterior and inferior wall of the left atrium. Rotors are established at the R1 and R2 zones. The high ApEn regions in Figure 8 are not specifically related to rotor activity. Some authors suggest that the standard ApEn parameters are not suitable for signals of fast dynamics, and that to solve these cases, the r and m parameters must be chosen from a larger set than the one proposed by Pincus [108, 114]. Following this idea, we designed an optimization process for the ApEn parameters obtaining the configuration ApEn(3, 0.38, 1000).Figure 8 (right) shows the corresponding ApEn in which the rotors R1 and R2 are highlighted by high ApEn values. Moreover, intermediate values of ApEn (green) are related to perturbations in conduction such as blockades, at the right inferior pulmonary vein, and shockwaves, at the zone below the coronary sinus. For additional details of this procedure and the results please refer to Ref. [115].

Figure 8.

ApEn electroanatomical maps for AF episode. The notations R1 and R2 correspond to the tip of rotors at the left atrial posterior wall and superior cava vein.

The fibrillatory activity presented in Figure 8 includes stable rotors that were characterized by the ApEn maps. We move now to the case in which the tip of the rotor meanders. To gain in temporal resolution, the parameter N was reduced to 500 points, and a nonoverlapping moving window was applied to each EGM of the model to obtain a time‐dependent ApEn value (Figure 9). We applied Pincus parameters ApEn(2,0.1,500) and optimized parameters ApEn(2,0.3,500). The analysis is performed over window of observation located at the left atrial posterior wall, in which a meandering rotor is generated. Figure 10 shows three consecutive frames of the episode (firs two columns) and the ApEn electroanatomical maps for Pincus (third column) and optimized (fourth column) parameters. For both parameter configuration, the high ApEn region changes as the tip of the rotor meanders through the observation window. The bottom row shows no rotor within the window that induces a reduction of ApEn for optimized parameters, while for Pincus parameters, high ApEn values remain.

Figure 9.

Nonoverlapping moving window procedure to obtain ApEn values varying in time.

Figure 10.

Three consecutive frames of propagation (first and second columns) of a meandering rotor and the Pincus parameters ApEn map (third column) and optimized parameters ApEn map (fourth column).

Under the assumptions of our computational model, we provide evidence that the hypotheses stated above: we are able to quantify fragmentation of EGM using the ApEn as a measure of regularity and to relate it with the tip of a stable and meandering rotors. Moreover, through an optimized version of ApEn parameters setup, other conduction anomalies can be identified. There are several works, with similar approach but with different tools of measurement of irregularity [54, 26, 88, 116, 50]. The tools used for the fragmentation analysis are mostly based on the calculation of the length of the cycle and the amplitude of the EGM, in correspondence with the definition of Nademanee [46]. Although the concept of CFAE established by Nademanee has been an important contribution to the study of AF, it may not describe the wide range of EGM fragmentations that occur in different cases. Therefore, we propose to extend the concept of fragmentation as a nonstatic and nonlinear phenomenon. The ApEn has already been applied in other studies for EGM analysis in AF [117], and in ventricular fibrillation [118] obtaining a single ApEn value for each EGM signal. We applied the Dynamic ApEn using a mobile window of 500 and 1000 points. An EGM of 4 s can provide up to 8 ApEn values, aimed to gain temporal resolution, and information about the behavior of EGM fragmentation.

The fragmented EGM were characterized by means of the ApEn, using the standard parameters suggested by Pincus and the parameters chosen from the proposed optimization method [115]. Both proposals reveal the relationship between CFAE and high values of ApEn, which is supported by Novak et al. [94]. Furthermore, EGM fragments with high ApEn values have been shown to be related to arrhythmogenic substrates, such as the rotor tip, blocking lines and the case of the coronary sinus area influenced by abrupt fiber direction, wave collision and passage from a narrow conducting zone to a wide but perpendicular zone. The relationship between CFAE and arrhythmogenic substrates has been reported in recent studies [61, 82, 87, 50, 49, 119–121]. However, how to differentiate the fragmented EGM according to the substrate that generates it? It has been observed that the ApEn maps, calculated from the standard parameters (ApEn (2, 0.1, 500) and ApEn (2, 0.1, 1000)), present fragmentation with high ApEn levels for the rotors R1 and R2, for the blockade at right inferior pulmonary vein and the coronary sinus region; however, they do not present significant numerical differences. This has also been observed by Navoret et al. [88], who establish a relationship between detected CFAE and the presence of wave and rotor shock, but they cannot differentiate them, using two‐dimensional computational models and cell cultures. They also report the development of an algorithm, which extracts five characteristics for the characterization of fragmented EGM. The algorithm is tested in two real databases: the first has EGM labeled as fragmented and nonfragmented. The algorithm shows good results in the classification. The second has EGM labeled as active fragmentation, whose ablation restored sinus rhythm, and fragmentation whose ablation did not restore sinus rhythm. The algorithm did not discriminate between both classifications [48]. However, the ApEn maps, calculated using the optimized parameters optimized (ApEn(3, 0.3, 500) and ApEn(3, 0.38, 1000)), assign values by ranges to areas of interest, the R1 and R2 rotors being the highest ApEn, followed by the intermediate values of ApEn in the zone of the blockade and the coronary sinus, and the smaller values of ApEn to the fibrillatory EGM of regular morphology. These results suggest that: the ApEn can solve the problem described by Navoret et al. [88]. On the other hand, if it is verified that the ablation guided by the ApEn maps restores the sinus rhythm, it solves the problem reported in Ref. [28]. Future work should focus on evaluating, in computational and experimental models, ablation guided by maps of dynamic ApEn.

Some authors have pointed out the disadvantages of ApEn: it is not a stable measure when the number of points N varies, and it has an inherent deviation from the real value, due to the inclusion of self‐comparison of segments [122, 123]. It has even been proposed a new statistical, the sample entropy (SampEn), as an enhancement of the ApEn [97]. However, it has been shown that both have a same behavior for time series of fast dynamics [110, 114, 124]. In our research, the possible instability that can present the ApEn due to the variation of points does not influence the results: specific families of ApEn have been defined, starting from the election of the r y m parameters, using time series of 500 and 1000 points. Although it has been observed that the range of ApEn values varies for each case, the information provided by the ApEn is embedded in the continuous scale that it results, given the degrees of irregularity that the fibrillatory EGM may present. This feature is important, since it offers adaptability to the wide range of EGM morphologies that may present different cases. In addition, it is not necessary to define intervals of fixed ApEn values for CFAE identification. It is only necessary to define that the ApEn scale, during an AF episode, is proportional to the presence of fragmentation, where the higher values of ApEn suggest the presence of rotors.

5.2. Multifractal analysis of EGM signals

Fluctuations in EGM signals are nonperiodic and exhibit nonlinear behavior. Fractal models could be used to describe this behavior through the identification of self‐similarity and scale invariance in the statistical properties of the signals. Fractal signals present self‐similarities and scale invariance properties that can be described by a single quantity, for example, the Hausdorff dimension or the Hurst exponent, and it is represented using a power law relationship. Fractal properties of physiological signals are not homogeneous, which means that local scaling properties change with time, and there is necessary to use different local Hurts exponent to describe the evolution of the system. Therefore, multifractal analysis could capture these changes in the global singularity distributions [101]. Then, the power law for multifractal behavior is written as follows:

$N_{\alpha }\sim {\mathrm{ϵ}}^{-f\left(\alpha \right)}$E4

where $N_{\alpha }$ corresponds to the number of balls with a singular exponent equal to some value of $\alpha$, necessary to cover a specific set, $\mathrm{ϵ}$ is the diameter of the balls that covered the set, and $f\left(\alpha \right)$ corresponds to the Hausdorff dimension. Note that in fractal analysis, the exponent is a scalar, while here the exponent is a function that contain different local Hurts exponent.

To compute the singular spectrum $f\left(\alpha \right)$, we used a method called multifractal detrended fluctuation analysis proposed and described in Ref. [125]. Using $f\left(\alpha \right)$, we calculate the h‐fluctuation index as is described by Orozco‐Duque et al. [101].

5.3. Location of simulated rotors using multifractal analysis

Multifractal analysis can be used to calculate features such as h‐fluctuation index and to build electroanatomical maps to located sites fractal or multifractal activity. This measure is related to the complexity of the EGM signals recorded over the atrial surface. In this case, we work with signals recorded from a simulated episode of AF in the three‐dimensional model described above. Figure 11A shows a color map built with the values obtained by the multifractal analysis according with Orozco‐Duque et al. [126]. Here, red dots are locating in regions where the rotor tip is presented in the model, and in the neighborhood of the rotor tip where the tip is meandering.

Figure 11.

(A) Three‐dimensional electroanatomical map using multifractal analysis. (B) Two‐dimensional electroanatomical map using ApEn. (C) Two‐dimensional electroanatomical map using multifractal analysis.

On the other hand, h‐fluctuation index was calculated in EGM signals recorded from the two‐dimensional simulated episode of fibrillatory conduction. For the sake of comparison, ApEn map was computed using the two‐dimensional model and the optimized parameters. These features were calculated for each EGM individually, and values were represented in a color map. Figure 11B shows the ApEn map and Figure 5 represented the MF map. ApEn map exhibits the behavior of the rotor tip and illustrates the dynamics in the vicinity of the pivot point. In the central point of the rotor, there is a blue point representing some organized signals; however, in the neighborhood of the tip, there are reds dots that represent signals with high fragmented activity. A start‐shaped pattern can be identified, which represents the movement of the singularity point.

In Figure 11C, multifractal map exhibits an interesting behavior because it captures the dynamic of the whole rotor in the atrial tissue, not only the neighborhood of the tip. One can note that the direction of the rotation is illustrated and can be interpreted using this map.

6. Combination of features in the same EA map

EGM signals exhibit different morphologies that have not been enough studied. An inadequate characterization of EGM morphologies has limited the success of EGM‐guide ablation strategies. Looking for a better description of different EGM patterns, some authors have proposed the combination of features to detect critical sites in AF. This approach has the advantage of use different information from the signals and detects patterns that could be associated with wave collisions, conduction block, or pivot points.

Ravelli et al. [127] used a logical digital map to combined two indexes, one based on the detection of activation rate, the cycle length; and another based on the analysis of similarities between activation waves, the similarity (S). Logical operator map can discriminate two morphologies: rapidly organized and rapidly fragmented. According with Kalifa et al. rapidly organized, EGM could be associated with localized source of AF, while rapidly fragmented could be associated with the neighborhood of a pivot point [52]. One limitation of the approach proposed by Ravelli et al. is that the computation of cycle length and S requires the segmentation of EGM signals and the detection of local activation waves (LAW). This process is not always feasible especially with high‐disorganized patterns.

Schilling et al. proposed an approach based on a classifier called Fuzzy decision tree to combine features and classify the EGM signals among four classes [128]. Classes are assigned according to a modified Well’s criteria, where class 0 corresponds to organized and nonfragmented signal, class 1 is assigned to signals with fragmented waves but periodic activity, class 2 corresponds to signals with fragmented waves with periodic and nonperiodic activity, and class 3 is assigned to signals with high frequency and continuous activity. A limitation for using supervised learning is that the classifier depends on the classes selected by the group who collected the training database. This could be critical and bias the results because there is not a complete understanding about EGM morphologies.

Unsupervised learning has been proposed to combine features and creates maps to show the distribution of EGM cluster according with conduction patterns. We have tested some clustering methods based on machine learning such as K‐means and Self‐Organized Maps (SOM) in previous work. To locate rotors over simulation in two‐dimensional models, the best performance was obtained using the combination of Shannon Entropy and the mean value of the EGM signal as an input of K‐means algorithm with three clusters or the SOM algorithm with four clusters. Figure 12B and C shows the result of K‐means and SOM applied to signals from the two‐dimensional model. The classifiers performance was calculated using the distance between the rotor tip location and the mean of the points grouped by the nearest cluster. These results evince the capability of the combination of simple features using unsupervised algorithms for rotor tip location.

Figure 12.

(A) Electroanatomical map using semisupervised clustering applied to signals from three‐dimensional model. (B) Two‐dimensional electroanatomical map using K‐means. (C) Electroanatomical map using SOM.

The issues related to the clustering approach are the selection of clusters number and the classification of clusters to identify the set of EGM related to a specific pattern. To overcome these limitations, Orozco‐Duque et al. [129] have proposed a semisupervised clustering approach to combine features. Semisupervised clustering (SC) is not limited to previously defined classes. The method selected was spectral clustering with an automatic detection of cluster number. In a previous work, it was tested the performance of SC to discriminated between four classes according with the scheme proposed by Schilling et al. [130]. We tested SC and its feasibility to located pivot points in simulated episodes of AF. We used unipolar signals acquired in an AF simulation using the three‐dimensional model of human atrial. Figure 12A shows a map built from the results of SC evaluation. Four clusters were detected; the cluster that represents regular signals is displayed in blue, and the cluster with the highest disorganized pattern is displayed in red. This cluster is located in the areas where the two rotors meander. The distribution of yellow and green clusters gives us an idea about cluster rotation.

7. Concluding remarks

Phase maps are, currently, the accepted tools for rotor characterization [131], by tracking the tip through the phase singularity. Although it is widely applied in computational studies [132, 133] and in vitro experimentation, such as optical mapping [134, 135], clinical application presents technical limitations: high spatial resolution and EGM require a preprocessing stage in which information of activation can be lost [136]. Here, we have presented how nonlinear features and clustering approaches provide information about the rotor dynamics during AF virtual episodes. Translated to a clinical context, these measures can be extracted from real EGM without needing special signal conditioning. However, our simulations provide high spatial resolution that places the same limitation as phase maps. In a recent report [137], we tested the influence of the spatial resolution over the ApEn electroanatomical maps in detecting rotors, using a two‐dimensional atrial fibrillation model. Our results indicate that the ApEn maps can identify the rotor tip with spatial resolutions close to those available in commercial mapping catheters. We also showed that a minor dependence of the ApEn maps on the virtual electrode array position, which implies that there is a transition of irregularity starting at the rotor tip and spreading to its surroundings. These findings encourage considering nonlinear features for EGM analysis. Although experimental validation is needed, further in silico studies are needed to enhance and characterize the behavior of these tools.

Nonlinear features such as ApEn and indexes calculated from multifractal analysis allow the construction of maps to display the distribution of EGM morphologies and to study the dynamic of fibrillatory conduction in the atrial surface. Additionally, application of clustering tools allows us to incorporate the information from different features within the same system for study the distribution of EGM clusters in the atrial surface. The use of unsupervised learning approach has the vantage that does not depend on a training specific dataset, which is an important feature considering the gaps in the knowledge about EGM morphologies. In AF simulated models, rotors were located by the proposed methodology; however, further observations and clinical studies are needed to associate marked sites with arrhythmogenic substrates in humans.