Atrial fibrosis identification with unipolar electrogram eigenvalue distribution analysis in multi-electrode arrays

Abstract Atrial fibrosis plays a key role in the initiation and progression of atrial fibrillation (AF). Atrial fibrosis is typically identified by a peak-to-peak amplitude of bipolar electrograms (b-EGMs) lower than 0.5 mV, which may be considered as ablation targets. Nevertheless, this approach disregards signal spatiotemporal information and b-EGM sensitivity to catheter orientation. To overcome these limitations, we propose the dominant-to-remaining eigenvalue dominance ratio (EIGDR) of unipolar electrograms (u-EGMs) within neighbor electrode cliques as a waveform dispersion measure, hypothesizing that it is correlated with the presence of fibrosis. A simulated 2D tissue with a fibrosis patch was used for validation. We computed EIGDR maps from both original and time-aligned u-EGMs, denoted as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}$$\end{document}R and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{R}^{\mathcal{A}}$$\end{document}RA, respectively, also mapping the gain in eigenvalue concentration obtained by the alignment, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \mathcal{R}^{\mathcal{A}}$$\end{document}ΔRA. The performance of each map in detecting fibrosis was evaluated in scenarios including noise and variable electrode-tissue distance. Best results were achieved by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{R}^{\mathcal{A}}$$\end{document}RA, reaching 94% detection accuracy, versus the 86% of b-EGMs voltage maps. The proposed strategy was also tested in real u-EGMs from fibrotic and non-fibrotic areas over 3D electroanatomical maps, supporting the ability of the EIGDRs as fibrosis markers, encouraging further studies to confirm their translation to clinical settings. Graphical Abstract Upper panels: map of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}^{\mathcal {A}}$$\end{document}RA from 3×3 cliques for Ψ= 0∘ and bipolar voltage map Vb-m, performed assuming a variable electrode-to-tissue distance and noisy u-EGMs (noise level σv = 46.4 μV ). Lower panels: detected fibrotic areas (brown), using the thresholds that maximize detection accuracy of each map


Introduction
Atrial fibrosis represents a structural anomaly of the atrial myocardium [1]. It is characterized by an altered extracellular matrix activity caused by fibroblasts [2], which alters the electrical conduction and excitability of the tissue. Fibroblasts activation and proliferation, as well as their secretion of extracellular matrix proteins, such as collagen, characterize fibrotic tissue [3]. This structural remodeling mainly occurs during the reparative process to replace damaged myocardial parenchyma [4]. In addition to this replacement process, others, as reactive fibrosis to a trigger as inflammation, have been recognized as responsible for fibrosis. This makes detection of fibrotic tissue even more difficult and suggests the need for more specific imaging tools and markers to detect and quantify fibrosis [5].
Atrial fibrosis has been observed to be correlated to atrial fibrillation (AF) [6]. Despite the fact that AF represents the most common cardiac arrhythmia, its trigger mechanisms are not yet fully understood [7] and its causal relationship with atrial fibrosis is still challenging [1]. On the one hand, fibrosis-induced remodeling creates a substrate promoting AF [1,6]; on the other hand, fibrosis can occur as a result of the electrical [1] as well as structural [4] atrial remodeling found in AF.
Atrial fibrosis is electrophysiologically characterized by low intracardiac electrograms (EGMs) amplitudes and conduction velocities [8], which may be measured with electroanatomical mapping (EAM) systems [9]. These allow displaying 3D voltage and activation time maps over a reconstruction of cardiac chambers anatomy and visualizing catheter position, so as to guide ablation procedures and treat arrhythmias with minimum radiation exposition [10]. Based on several studies, consensus exists on the choice of 0.5 mV as the threshold value of bipolar EGMs (b-EGMs) peak-topeak amplitude to discriminate fibrotic areas in the atrium during sinus rhythm [11]. However, this procedure brings along several limitations. First of all, a peak-to-peak voltage measure disregards morphological and temporal information contained in the signal. Nevertheless, this reflects the possible presence of underlying abnormalities in the atria [12]. Secondly, bipolar voltage mapping may be influenced by other technical factors not related to the substrate, including the relative orientation between the recording electrode pair and wavefront propagation direction, electrode size, interelectrode distance and b-EGMs filtering [12]. Bipolar voltage is also affected by the tissue-electrode contact, whose maintenance may be challenging in anatomically difficult sites (e.g., the pulmonary veins). Third, the definition of low-voltage areas has not been subjected to a standardization procedure and the voltage threshold has never received an histological validation [12].
In recent years, more attention has been paid to the role of fibrotic tissue on the initiation and perpetuation of AF than on its effects on the morphology of the EGMs [13]. In this sense, not many intracardiac signal processing methods based on EGM features have been proposed to detect fibrosis. Some studies have revealed the relationship of EGM morphology and tissue alterations, including ablation lesion formation [14]. Others have introduced a method to characterize the different fibrotic textures based on EGM fractionation due to the incidence of wavefront direction [15]. All these works have used in silico experiments for validation.
In this paper, we hypothesize that the waveform dispersion of neighbor unipolar EGMs (u-EGMs) is correlated with the presence of atrial fibrosis. Therefore, the aim of the work is to propose eigenvalue-based indices of waveform dispersion to identify fibrotic areas. They take into account the spatiotemporal relations of the u-EGM waveforms and overcome the limitations of the use of b-EGM voltage thresholding for the detection of fibrosis. Resulting maps, called eigenvalue dominance ratio (EIGDR) maps, were computed with a simulation setup, using two sizes of nearby electrode arrangements, referred as cliques (2× 2 and 3 × 3 electrodes), and three catheter orientations with respect to the tissue preferential direction (0 • , 30 • and 45 • ). The ability of each map was evaluated in detecting a fibrosis patch included in the simulated tissue. As a proof of concept, EIGDR values were also computed in clinical u-EGMs from patients, in four-and five-electrode cliques. They were correlated with the presence/absence of fibrosis based on late gadolinium enhancement-magnetic resonance imaging.
This paper is organized as follows: Section 2 presents the methodology and the datasets used for its validation. Section 3 contains the obtained results, whereas Sections 4, 5 and 6 include the discussion, limitations (with references to challenges and future perspectives) and conclusions of the work, respectively.

Atrial model
We simulated an atrial tissue slice of 4 × 4 cm by dividing it into adjacent square elements whose centers were separated 0.1 mm. Within the tissue slice, a circular patch with a diameter of 2 cm was defined, whose center was at the center of the 2D tissue, as shown in Fig. 1(a). Inside the patch, a fully transmural (from endo-to epicardial layer) pattern of diffuse fibrosis was randomly defined following a uniform distribution. The Maleckar model for myofibroblasts [16] was assigned to 20% of the nodes within the circular area. Although the percentage of atrial fibrosis is very patientdependent, the density of 20% represents the threshold value between stage II and stage III according to the Utah classification [17] and therefore considered a realistic percentage for this study.
The cell model assigned to all the non-fibrotic nodes was a variant of the Courtemanche myocyte model [18] accounting for the atrial electrophysiological characteristics experimentally observed in the left atrium (LA) and the persistent AF (cAF) induced remodeling [19]. This electrical remodeling was introduced through the variation of the maximum conductances of the transient outward potassium current (I to ) , the L-Type calcium current (I CaL ) , the inward rectifier potassium current (I K1 ) , the ultrarapid outward potassium current (I Kur ) and the slow delayed rectifier potassium current (I Ks ) (see Table 1), as in previous computational studies [19,20]. Additionally, the diffusion tensor was adjusted to reproduce the conduction velocity in the LA, and was reduced by 30% in all elements of the patch with at least one fibroblast node [21], similarly to what was performed in [22].
Simulations were run in ELVIRA software [23]. The monodomain formulation was solved using the operator splitting numerical scheme with a constant time step t = 0.01 ms and a spatial resolution x = 0.1 mm. Recording electrodes were distributed over the simulated anatomy mimicking a L × L high-density multi-electrode array (MEA) (L = 15), where the inter-electrode distance was d = 2 mm. The simulated electrode grid was located so that its central electrode corresponds to the center of both the tissue slice and the fibrotic patch, and was rotated by an angle Ψ , Ψ ∈ {0 • , 30 • , 45 • } , with respect to the tissue fiber direction. Fig. 1(b) shows the three MEA orientations over the simulated tissue geometry. Two generic cliques of different sizes were considered in this work as depicted in Fig. 2(a) and (b), forming four and nine-electrode square arrangements, respectively. Each clique location is referred to the lower left electrode, indexed as (i, j), within the MEA. The rest of electrodes in the clique are numbered from left to right and bottom to top, thus corresponding to locations (i , ⋯ , 13} ) for the 3 × 3 arrangements, respectively.

Synthetic signals
Unipolar EGMs, u i,j (n) , were calculated as originated by the passage of the propagation wavefront by electrodes located at sites (i, j) of the MEA. We computed them in a volumetric tissue-blood model, assuming a temporal resolution of 1 ms, as done in [30]. First, extracellular potentials were obtained by an approximation of the bidomain formulation, considering the tissue immersed in a non-conductive bath. In order to reproduce fibrosis effects, inside the simulated fibrotic patch the single cardiomyocyte under cAF conditions was coupled with a randomly variable number of fibroblasts within the patch. Fig. 3(a) shows the action potentials registered in two different cardiomyocytes from the mesh, one outside the fibrotic patch and the other inside the fibrotic patch, where it is coupled with two fibroblasts. Electrical remodeling induced by cAF produces a 55% reduction in duration measured at 90% repolarization (248 vs. 111 ms), in concordance with experimental data [29]. Fibroblast coupling with cAF cardiomyocytes makes resting potential less negative (83 vs. 78 mV) and elongates the duration measured at 90% repolarization (111 vs. 120 ms). Then, u-EGMs were solved in the entire domain by the governing equations for a solid volume conductor and its boundary conditions at the tissue-blood interface.
In order to take into account a realistic scenario, where the distance from the tissue may not be constant through the electrodes in the clique, we considered that the electrode-tissue distance, i,j , varies following a normal distribution with mean = 1 mm and standard deviation = 0.1 mm. Two thousand different random configurations were simulated, where the distance of each electrode to the tissue was randomly and independently chosen following that distribution.
Synthetic u-EGMs were computed with sampling frequency of 1 kHz, duration of 0.5 s ( N = 500 samples), and including one single activation (depolarization plus repolarization).
Simulated signals were corrupted with noise excerpts obtained from real u-EGMs, as previously done in [8]. Two thousand different noise segments were extracted from u-EGMs recorded with a multi-electrode PentaRay ® catheter (Biosense-Webster, Inc., Diamond Bar, CA, USA) at intervals with no recorded EGMs. All noise segments were normalized to have standard deviations v ∈ {0.0, 5.8, 11.6, 23.

Clinical data
Intracavitary u-EGMs recorded during sinus rhythm with a PentaRay ® catheter (Biosense-Webster, Inc., Diamond Bar, CA, USA) were used to evaluate performance of EIGDR approach with real signals. Clinical data were obtained from a patient with cAF, a slightly dilated LA (26 cm 2 ), a left ventricular ejection fraction of 58% and treated with anticoagulation and flecainide, registered at the Hospital Clínic, Barcelona, Spain. The data acquisition protocol was reviewed and approved by the Hospital Clinic Ethical Committee (Ethics approval number: HCB/2019/0881). The patient was informed and signed the consent form. A total of 758 mapping points (or catheter sites) were acquired at the anterior, posterior, lateral and septal wall,  [24] as well as the left atrial appendage and the pulmonary veins of the LA (left atrial regions) using the CARTO ® 3 EAM system (Biosense-Webster, Inc., Diamond Bar, CA, USA), so as to reconstruct a real-time 3D anatomical map before the ablation procedure. Signals were acquired from the 20 poles distributed among the five branches of the catheter, Fig. 4, characterized by consecutive inter-electrode spacings d, at each branch, of 2, 6, and 2 mm, resulting in 20 u-EGMs associated at each mapping point. In order to evaluate performance of EIGDR-based markers, 38 catheter positions were selected by an operator-dependent visual approach. Specifically, we visually identified and manually annotated nineteen points clearly assignable to fibrotic and another nineteen points to non-fibrotic areas on the anatomic map. In order to guide this decision, the corresponding magnetic resonance image (MRI) was used as reference.
Three of the catheter positions selected at fibrosis showed poor electrode-tissue contact for one or more splines of the catheter. In addition, mapping points located at borderline areas between fibrotic and healthy tissues over the MRI and/ or the bipolar voltage map, which showed bipolar amplitudes unhealthy (< 0.1 mV) or poor healthy (< 0.5 mV), have been excluded from the analysis.
EAM data and MRI were co-registered with the ADAS 3D Medical imaging software (ADAS-3D, Barcelona, Spain), as shown in Fig. 5. That co-register was performed by manually selecting some landmarks (between six and ten) in specific areas (such as the pulmonary veins and the atrial appendage) of the meshes. In order to determine how pathological the tissue is, the methodology described in [32] was used. Following the image intensity ratio (IIR) based thresholding, a color-coded 3D mesh was automatically generated, showing in blue the healthy tissue (IIR<1.2) and in red the dense fibrosis (IIR>1.32) (Fig. 5).

Eigenvalue analysis
Unipolar signals in the clique can be compactly represented by the following N × K matrix, K ∈ {4, 9} , as in [33]: where the k-th column, k , contains the samples of the unipolar signal u k (n) , modeled later in Section 2.6: We propose and assess EIGDR values from the clique u-EGMs in (1) to detect fibrosis. The N eigenvalues were obtained from the following correlation matrix estimate within each clique: The matrix u is the intra-signal sample correlation matrix, whose eigenvalues reflect the degree of morphological variability among the signals in the clique.
For each clique, the ratio R of the largest (i.e. dominant) eigenvalue 1 and the remaining ones was estimated: By comparing the first eigenvalue to the sum of the others, we are able to quantify the u-EGM energy percentage which can be explained by the shape of the first eigenvector. The ratio R would be much higher than one if all u-EGMs were essentially identical to each other and all the waveforms can be explained with just the shape of the first eigenvector (i.e., low morphological variability). On the contrary, when waveform dispersion appears, the eigenvalues from 2 to N become higher, thus reducing the ratio R.

Wave alignment
Eigenvalues of ̂ u were computed from the original u-EGMs and after intra-clique time alignment, proposed to compensate the effect that different activation wavefront arrival times have on the EIGDR. In this case, unipolar signals were aligned according to Woody's iterative procedure [34,35]: 1) at first l-th iteration, l = 0 , the relative time delay ̂k ,0 between each unipolar signal u k (n) and the u-EGM with the highest peak-to-peak amplitude within the clique u max (n) (assumed as initial reference signal) was estimated by maximizing their cross-correlation: 2) from the relative time delays, the average of the shifted signals u k (n −̂k ,0 ) within the clique was calculated: 3) in each l-th iteration, l > 0 , the cross-correlation between each u k (n) and the average signal obtained in the previous iteration ū l−1 (n) (assumed as updated reference signal) is maximized to find the updated relative time delays ̂k ,l : This process is repeated iteratively until the delay estimates no longer change. Eigenvalues of the intra-signal correlation matrix of the aligned u-EGMs were calculated in the same way as for their non-aligned version, thus leading to the formulation of the ratio R A (where the upper index denotes that the ratio comes from aligned u-EGMs within the clique) analogous to (4).

Unipolar signal modeling
The u-EGM signals within a clique located at (i, j) position within the MEA are indexed as u k (n) , k ∈ {1, ⋯ , K} , K ∈ {4, 9} . Assuming a plane wave propagation, the different electrodes in a clique are activated at different times, and therefore u-EGMs in the clique u k (n) will be delayed versions among them, plus noise and non-homogeneous components. Therefore, they can be modeled, analogously to [36] for misaligned signal ensembles, as: where: • s(n) is the u-EGM activation signal component assumed to be space invariant in the case of a plane wave propagation and homogeneous tissue free of fibrosis. Its energy is denoted as E s . • k is the delay of the k-th u-EGM s(n − k ) with respect to a time reference within the clique, introduced later. Delays k are zero-mean and characterized by their variance in the normal tissue 2 . In fibrotic areas, the reduction in conduction velocity with respect to healthy tissue increases the variance of the k up to 2 2 , where the factor > 1 in fibrosis and = 1 in non-fibrotic tissue. • k is a parameter accounting for u-EGM amplitude reduction between fibrotic, k < 1 , and healthy, k = 1 , tissues. It is modeled as a random variable with mean and variance 2 ( = 1 , 2 = 0 in healthy tissue; < 1 , 2 > 0 in fibrosis). • f k (n) is a zero-mean fibrotic signal component across the clique with variance 2 f . In healthy tissue f k (n) = 0. • v k (n) is the zero-mean noise component at the k-th u-EGM, with variance 2 v , Gaussian, white and uncorrelated with k and f k (n).
Four different scenarios for each u k (n) were considered in this study, as already proposed in [37], with/without fibrosis and with/without prior alignment. Their approximate theoretical eigenvalues and EIGDR were derived following parallel methodology to that used in [36] for repetitive signal ensemble alignment, as detailed below.

Prior alignment with no fibrosis
For perfectly aligned signals without fibrosis, i.e., u k (n) = s(n) + v k (n) , the intra-signal correlation matrix is given by: where is the N × N identity matrix and data vector , T , is easily shown to be proportional to the first eigenvector of u , whereas the remaining eigenvectors are chosen arbitrarily as long as they are orthogonal to the first. The eigenvalues are given by: where E s = T is the signal energy. For real u-EGM, signal energy is much larger than noise energy, E s >> N 2 v , and N >> 1 , resulting in an EIGDR of:

No prior alignment and no fibrosis
We now analyze the case of raw u-EGMs in the clique where misalignment of s(n) is assumed to be present in each kth u-EGM, u k (n) = s(n − k ) + v k (n) . Then, to estimate the eigenvalues we can approximate u k (n) , for small k , as [36]: where s � (n) and s �� (n) denote the first and second derivative of s(n), respectively. The intra-signal correlation matrix can be expressed as: where ′ and ′′ are the vector counterparts of s � (n) and s �� (n) , respectively. It can be shown that the eigenvalues of u are approximated by [36]: where E s � = �T � is the derivative signal energy. The resulting EIGDR is approximated by: Note that when = 0 (i.e., perfect alignment) this equation becomes equal to (11).

Prior alignment and fibrosis
When u-EGMs are first aligned in fibrosis zones, each of them can be modeled as u k (n) = k s(n) + f k (n) + v k (n) . The correlation matrix results in: and their corresponding eigenvalues are: which lead to: where the parameters used have already been introduced at the beginning of this section and upper A and lower F indices in R A F denote that the calculations are obtained from intra-clique aligned u-EGMs in fibrotic tissue, respectively.

No prior alignment with fibrosis
When raw u-EGMs come from fibrotic areas, each of them can be modeled as u k (n) = k s(n − k ) + f k (n) + v k (n) . In this case, the delay k has larger standard deviation than in non-fibrotic areas, with the extra delay controlled by a multiplicative factor , thus resulting in the following correlation matrix: which is similar to the case without fibrosis in (13) but with different proportionality factor and two different random components, corresponding to noise and fibrosis. The eigenvalues will similarly be approximated by: The corresponding EIGDR is approximated by: A summary of eigenvalues and the EIGDR are reported in Table 2 for the four scenarios.
Note that if the inter-signal correlation matrix, been computed rather than the intra-signal correlation matrix u , for the more general case with fibrosis and misalignment, and following a derivation parallel to the one presented in [36], the eigenvalues would have resulted in: When computing the EIGDRs for this matrix the results are approximately the same as for u . In practice, the matrix is estimated as: rather than with the theoretical expectations. From these estimates we observe that matrix ̂ • u is full rank while matrix ̂ u (estimated as in (3)) is not ( N > K ), circumstance that does not represent a limitation since no matrix inversions are required. In addition, as shown in [33], for the dataestimated autocorrelation matrices with , again showing equivalence of the EIGDR ratios for both matrices when estimated from the available data. Therefore just computational considerations can advise to use one or the other.

EIGDR-based fibrosis markers
According to the previous model, three main differential effects on signal shape, amplitude and arrival times are expected to occur in fibrotic as compared to non-fibrotic tissue: 1) higher morphology dispersion represented by the u-EGM signal component f k (n) and quantified by 2 f ; 2) lower and less homogeneous signal amplitudes represented by k and quantified by < 1 and 2 (( 2 + 2 ) < 1); 3) larger inter-signal misalignment within the clique as a result of slowed conduction and represented by delays k with enlarged variance ( 2 2 , > 1 ) relative to healthy tissue ( 2 ).
In order to determine which ratio to use for discriminating F and NF tissues, we first compare the EIGDR computed with no prior alignment of u-EGMs and no fibrosis, R , with the case with fibrotic tissue, R F , see Table 2. We observe that R F < R , since the numerator in R F is smaller than in R as a result of being larger than one, while the terms in denominator are larger, 2 > 1 , 2 f > 0 and ( 2 + 2 ) < 1 , as a combination of the three effects introduced by fibrosis. This result suggests the use of the ratio R , which becomes R F in fibrosis, with a thresholding strategy to discriminate if the clique is at fibrotic or healthy tissue. The same analysis can be done when comparing the EIGDR with prior alignment of u-EGMs and fibrosis, R A F , with respect to its counterpart with no fibrosis, R A . In this case, only the terms 2 f > 0 and ( 2 + 2 ) < 1 are responsible of the difference, since 2 has already been compensated for with alignment, resulting in R A F < R A . This suggests that R A may also be used as a thresholding strategy to discriminate fibrotic from healthy tissue.
In order to study which of the two options, R or R A , is more sensitive to fibrotic tissue characteristics, we analyze how the difference between R and R F evolves by varying 2 . For that purpose, we compute the ratio ΔR F , which under the proposed signal model can be approximated by: Its partial derivative with respect to 2 is: Typically 1 − 2 2 + 2 > 0 in fibrosis since conduction velocity reduction is less prominent ( ≈ 2 for high fibrosis, [38,39]) than voltage attenuation ( ≈ 0.3 [12]) and consequently the product 2 2 < 1 . This results in meaning that the lower the misalignment 2 the larger ΔR F , implying higher EIGDR differences between fibrotic, R F , and healthy, R , tissue. This justifies the advantage of prealigning u-EGMs in the cliques before EIGDR calculations, since the higher the misalignment 2 , the lower ΔR F and consequently the capacity of R to discriminate between fibrosis and non-fibrosis, and suggests that R A is better suited fibrosis marker than R. Table 2 Signal models for non-aligned (NA) and aligned (A) u-EGMs at non-fibrotic (NF) and fibrotic (F) areas, with their respective eigenvalues k and eigenvalue dominance ratios EIGDR u-EGM model Alternatively, we considered the ratio ΔR A of EIGDR computed from perfectly aligned u-EGMs with respect to misaligned original u-EGM signals, representing the gain in eigenvalue concentration produced by alignment: This expression has been estimated for the more general case including fibrosis ( 2 f > 0 ), so expressions R A F and R F are used. Nevertheless, it can certainly be computed at cliques on any tissue, fibrotic or non-fibrotic, and its theoretical value when no fibrosis is present can be retrieved from (26) just by making f = 0 , =1 , = 1 , and = 0 . Sensitivity of ΔR A to fibrosis has been analyzed by deriving (26) with respect to the fibrosis-induced parameters. Therefore, deriving with respect to 2 f to see how ΔR A depends on the level of fibrosis, we obtain: f increases and suggesting the possibility of using ΔR A as a fibrosis marker, like R and R A . This behavior is also corroborated by computing the derivative of ΔR A with respect to 2 , taking into account u-EGM amplitude reduction in fibrosis: This expression results > 0 for small delays k , confirming that the larger the fibrosis (i.e., the smaller 2 ), the smaller Nevertheless, under the same assumptions, the derivative of ΔR A with respect to 2 results in ΔR A 2 > 0: meaning that the larger the reduction of velocity due to fibrosis (i.e., the higher ), the greater ΔR A , thus showing an opposite effect.
However, as already said before, we expect that fibrosis effects on u-EGM amplitude and morphology, expressed by 2 and 2 f , respectively, are much more marked than those on (26) , conduction velocity given by 2 [40]. Therefore, we expect the index ΔR A to be reduced when fibrosis is more severe. According to this analysis, three different EIGDR-based metrics revealed to be sensitive to the presence of fibrosis and therefore suitable to distinguish between fibrotic and nonfibrotic areas: R , R A , which can be interpreted as measurements of the shape homogeneity of the u-EGMs before and after time alignment, respectively, and the ratio between both, ΔR A . Maps of R , R A and ΔR A were computed by processing the complete MEA in the two clique sizes considered in the simulation study. Each map consists of color-coded pixels, representing EIGDR value at each clique. The 2 × 2 configuration provides one EIGDR value for each square group of four electrodes with diagonal vertices at (i, j), and (i + 1, j + 1), i, j ∈ {1, ⋯ , 14} , giving a total of 14 × 14 pixel maps for each proposed marker. The 3 × 3 clique provides one EIGDR value at each squared group of nine electrodes with diagonal vertices at (i, j) and (i + 2, j + 2) , i, j ∈ {1, ⋯ , 13} , resulting in maps of 13 × 13 pixels, for R , R A and ΔR A .

EIGDR with variable electrode-to-tissue distance
When we introduce variable electrode-to-tissue distance, we need to modify the model by replacing s(n) with k s(n) , being k a random variable with mean = 1 and variance 2 indexing all the i,j within the clique. Similar analysis as in previous subsections leads to obtain: which just introduces a multiplying factor, (1 + 2 ) , with respect to the ratios in (11) and (18), equal in both ratios, so preserving the fibrosis stratification value of R A biomarker. This occurs in contrast to b-EGM peak-to-peak marker, V b i,j , which is largely modified by the variable electrode-to-tissue distance, but a random way at each electrode, reducing its value as a stratification marker. Analogously: which approximately result in the same multiplying factor (1 + 2 ) with respect to ratios (15) and (21). Note that for small , N 2 v >> 2 E s ′ and then the approximation of a multiplying factor relating fixed with variable electrode-to-tissue distance holds. This again shows that R preserves the fibrotic stratification value in variable electrode-to-tissue distance situations. The same analysis also applies to the ratio ΔR A .
Also note that in presence of more complex fibrillatory propagation patterns, changes occurring in the u-EGM morphology from electrode to electrode can initially be thought as a planar wave propagating in different directions. This will also result in an extra k-dependent amplitude component into the s k (n − k ) signal in the model of (8), depending of the angle of the planar wave, and thus also evidence not to largely affect the EIGDR.

EIGDR in real data from PentaRay®
Values of R , R A and ΔR A were also computed within the four-and five-electrode cliques considered at each mapping point acquired by the PentaRay ® catheter. In order to quantify the atrial activity related dispersion, an atrial depolarization window of 100 ms fixed length ( N = 100 ) was extracted from the last recorded activation at each recording site. Therefore, the proposed EIGDR-based markers were calculated using windowed signals in each clique, aligning them when required as explained in Section 2.5.

Voltage-based fibrosis markers
We also considered bipolar voltage maps based on the peak-topeak amplitudes , j ∈ {1, ⋯ , 14} . These peak-to-peak amplitude-based maps were considered and their performance for fibrosis detection were compared against EIGDR maps. Each color-coded pixel bipolar map presents the same resolution as 2 × 2 cliques EIGDR maps, providing 14 × 14 pixels when processing the whole MEA.
Regarding clinical data, for each mapping point we derived b-EGMs along the PentaRay ® catheter branches from filtered u-EGMs. Peak-to-peak amplitudes were computed using atrial depolarization windows extracted from the last recorded activation of b-EGMs at each bipole.

Performance assessment for fibrosis detection
Both EIGDR and bipolar mapping strategies were estimated and tested for each noisy realization u q i,j (n) , i, j ∈ {1, ⋯ , 15} , q ∈ {1, ⋯ , 2000} considered in this study. For each map type, results are reported by aggregating the three different MEA orientations. This aggregated version represents a scenario where the relative angle of the propagation direction with respect to the catheter was not known a priori, thus being more realistic.
In order to quantitatively evaluate the ability of the different maps as markers for fibrosis detection, i.e. in discriminating pixels associated to the fibrotic patch from those related to non-fibrotic tissue, receiver operating characteristic (ROC) curves were used. Two ground-truth masks were created for that purpose, with the resolution of 14 × 14 and 13 × 13 maps, by labeling whether a clique (or an electrode pair in case of bipolar maps) lies within the fibrotic or the non-fibrotic area. In a first analysis, the 14 × 14 ground-truth mask was created by assigning value 1 if the four electrodes within a 2 × 2 clique lie in the fibrotic area, and value 0 if the four electrodes lie in the non-fibrotic area. In a similar way, the 13 × 13 groundtruth mask was created by considering if the nine electrodes within a 3 × 3 clique fully lie or not in fibrotic/non-fibrotic tissue. Cliques with some electrodes inside and some outside the fibrotic patch were not labeled and therefore discarded in the evaluation. The two ground-truth masks used in this study are shown in Fig. 7(a) and (b), for evaluating 14× 14 maps (both EIGDR and bipolar) and 13× 13 maps, respectively. Then, in a further analysis, two binary 14 × 14 and 13 × 13 ground-truth masks including those mixed cliques with electrodes inside and outside fibrosis region were considered for the evaluation. For that purpose, cliques whose distance between their central point and the center of the fibrotic patch was shorter than the radius of the patch were labeled as fibrotic. On the contrary, when this distance was longer than the radius, corresponding cliques were classified as non-fibrotic. For each EIGDR and bipolar mapping strategy, ROC curves were computed by varying the threshold for fibrosis identification, obtaining sensitivity and specificity in the detection of the fibrotic area [41]. In this context, true positive denotes the number of cliques correctly identified as fibrotic, false negative represents the number of missed cliques in the fibrotic area, true negative stands for the number of cliques correctly identified as non-fibrotic and false positive is the number of cliques incorrectly detected as fibrotic. The maximum detection accuracy (ACC ), defined as the highest number of correctly identified cliques (fibrosis or non-fibrosis) divided by the total number of assessed cliques, was used as a measure of the overall fibrosis detection ability of each map. Values of ACC , as well as of the threshold corresponding to ACC , were computed for each map aggregated version considered in this work. Averaged results over the noisy realizations were then computed and evaluated for performance measurements.
In the clinical data analysis, at each F and NF mapping point, median values of the three ratios R , R A and ΔR A were calculated over the six cliques considered. Analogously, the median and maximum values, V b and V b-m , respectively, among the five peak-to-peak bipolar amplitudes computed along the catheter branches and associated with the bipoles closest to the its center ((3,4), (7,8), (11,12), (15,16) and (19,20), according to Fig. 4), were computed at each mapping point in fibrosis and healthy tissue. In order to compare markers between F and NF tissues, median and interquartile range (IQR) over all the EIGDR-based indices and bipolar amplitudes were computed at both F and NF points, separately, as well as the p-values of the right-tailed Wilcoxon rank-sum test referring to the comparison of the metrics between the two areas. Finally, median and IQR of EIGDR indices and bipolar amplitudes were calculated over the six cliques and the five innermost bipoles, respectively, of all mapping points considered, at both F and NF areas, as well as their p-values referring to the global comparison of the metrics between those F and NF areas.

Analysis of simulated data
An example of the mapping strategies (computed with 3 × 3 cliques) for Ψ = 45 • , variable electrode-to-tissue distance and without noise, is shown for EIGDR and bipolar maps at upper panels in Fig. 8(a) and (b), respectively. When noise is added at level of v = 46.4 V , results for one of the two thousand noisy realizations are presented at Fig. 8(c) and (d). In the lower panels, the fibrotic areas identified by using the thresholds that maximize the ACC are shown for each mapping strategy. Blue (brown) color inside the circle encompassing fibrotic patch denotes false negative (true positive), while outside denotes true negative (false positive) detection, respectively.
We reported values of R , R A and ΔR A computed from noisy ( v = 46.4 V ) u-EGMs in the non-fibrotic clique Results in this example illustrate that EIGDR maps performed from noise-free time-aligned u-EGMs, R A and ΔR A , plotted at central and rightmost columns in Fig. 8(a), respectively, present fibrosis detection performance comparable to bipolar maps obtained as the maximum voltage of both MEA directions, V b-m , shown at rightmost column in Fig. 8(b). However, when u-EGMs are affected by noise, EIGDR maps (upper row at Fig. 8(c)) are more robust than bipolar maps (upper row at Fig. 8(d)), being R A the one showing the best fibrosis discrimination performance.
ACC values of all mapping strategies considered in this study are summarized in Table 3, assuming fixed or variable distance between MEA and tissue, and five different noise levels (reported as standard deviations v and average peakto-peak amplitudes V pp,v ). For each map, thresholds having these maximum detection accuracy values were also reported in Table 4, where they were given in voltage units in case of b-EGM amplitude-based maps. Both ACC and threshold values were calculated and reported by aggregating the three catheter orientations with respect to the propagation direction. Despite this, bipolar voltage maps exhibit performance

Fibrosis
Fixed electrode-to-tissue distance (FD) Variable electrode-to-tissue distance (VD) strongly dependent on the relative orientation between MEA and propagation direction (e.g., ACC = 68.7% and ACC = 90.8% for V b-x and V b-y , respectively, for fixed catheter-to-tissue distance and v = 0.0 V ). For v = 46.4 V , V b-y and V b-m are more affected by noise than EIGDR maps, both for fixed and variable electrode-to-tissue distance. The selection of the thresholds for the EIGDR implies another challenge. This can be addressed by observing that values reached by R A at healthy tissue (11) can be obtained as the ratio between the estimated energy Ê s and N times the estimate of the noise variance ̂2 v , R A =Ê s ∕(N̂2 v ) . The value of Ê s can be estimated from the data (e.g., by averaging the energy of the u-EGMs in the clique), while the noise variance ̂2 v can be estimated as the u-EGM signal variance in areas electrically silent. The threshold, T , can be fixed to a value T =Ê s ∕(N̂2 v ) − Δ , where Δ will control the trade-off between required sensitivity and specificity: small Δ will provide high sensitivity and low specificity, and the reverse for large Δ . Note that from results reported in this subsection for a non-fibrotic site such as (i, j) = (3, 3) with noise v = 46.4 V , we measured R A = 7.76 , while the optimum threshold reported in Table 4 for this noise level is 3.5, corresponding to a Δ = 4.26 , which can be taken as a reference value. Similar analysis for the same noise level gives Δ values estimates of 1.7 for R and 0.39 for ΔR A , as quantities to subtract to the estimates of R and Δ R A at non-fibrotic areas to derive usable threshold values in real clinical settings. These estimates will require additionally an estimate of the ̂ and Ê s ′ , see Eqs. (15) and (26), which can be computed, e.g., from the standard deviation of estimated delays in a clique and from the derivative of the aligned and averaged u-EGMs in the clique, respectively. Bipolar voltage map V b-m identifies the fibrotic area with an ACC of 96.2% when distance between MEA and tissue is fixed and u-EGMs are not affected by noise ( v = 0.0 V ). Nevertheless, this performance reduces to ACC = 92.5 ± 1.1% when the electrode-to-tissue distance is variable, and further when increasing noise level, reaching values 86.9±1.1% and 86.1±1.2% for the highest noise level ( v = 46.4 V ), in case of fixed and variable electrode-to-tissue distance, respectively. On the other hand, R A performed with 3 × 3 cliques is more robust to the effect of variable distance than V b-m , presenting ACC = 92.1% and ACC = 92.3 ± 0.8% for fixed and variable catheter-to-tissue distance, respectively. In addition, it achieves ACC = 94.2 ± 1.6% when v = 46.4 V , for both fixed and variable distance scenarios, being consistent with example in Fig. 8. The same behavior has been observed when studying the three MEA orientations separately. For the highest noise level under test, and with 3 × 3 cliques, R A achieves greater ACC values than V b-m . In particular, R A reaches 95±2%, 95±3%, and 95±3% for Ψ = 0 • , 30 • , and 45 • , respectively, both with fixed and variable electrode-to-tissue distances, while V b-m reaches 88±2%, 89±2%, and 90±2% with fixed distance, and 87±2%, 88±2% and 89 ±2%, with variable distance, for Ψ = 0 • , 30 • , Table 4 Thresholds corresponding to the ACC values reported in Table 3, presented as mean ± standard deviation except for fixed electrode-to-noise distance and v Variable electrode-to-tissue distance (VD) and 45 • , respectively. If the evaluation is performed without exclusion of cliques which have electrodes inside and outside the fibrotic patch, R A still provides higher ACC (83.0 ± 1.5%) than V b-m (81.2 ± 1.21%), in the largest noise contamination and with variable electrode-to-tissue distance. Moreover, when u-EGMs are not affected by noise, ACC goes from 80.2 % to 80.1 ± 0.6 % for R A while ACC reduces from 90.8 % to 87.2 ± 1.0 % for V b-m , from fixed to variable electrode-to-tissue distance. Table 5 shows the median values and IQR of the different biomarkers computed at nineteen mapping points at fibrotic (F) and other  [3.12] for R A and ΔR A at non-fibrotic and fibrotic tissue, respectively. These results are consistent with the theoretical model and overtake fibrosis discrimination performance of V b (Wilcoxon rank-sum test, p-value<0.05).

Clinical significance of the work
Detection of atrial fibrosis is capital for guiding catheter ablation strategies in AF. The typical intra-procedural assessment of atrial fibrosis by means bipolar voltage thresholding presents well-established limitations related to catheter-wavefront orientation, catheter-tissue contact, electrodes size and interelectrode spacing, thus limiting its reliability as surrogate of fibrosis. Besides this, it is well-known that when using threshold-based approach, EGM morphology information and time relationship among adjacent electrodes are missing. Despite late gadolinium enhancement-magnetic resonance imaging represents the only non-invasive tool for atrial fibrosis diagnosis, its reproducibility remains under debate [42], as well as its utility in clinical settings [43].
In this work, we proposed u-EGMs eigenvalue dominance ratios (EIGDR) to quantify voltage waveform dispersion and investigated their performance as markers in discriminating fibrotic and nonfibrotic areas, by using a 2D simulated tissue including diffuse fibrosis. The hypothesis behind this approach is that underlying fibrosis in the atrium is reflected not only in the reduction of the waveform amplitude but also in the increased inter-signal dispersion in cliques of nearby electrodes, and that this dispersion will be insensitive to electrode-to-tissue distance as opposite to b-EGM amplitudes.

Performance evaluation of fibrosis markers with simulated data
We analyzed maps computed from noise-free u-EGMs, as well as from u-EGMs corrupted by homogeneous noise  levels. As a first step, the distance between each electrode of a square MEA and tissue was assumed to be fixed at 1 mm. Then, in a further analysis, that distance was assumed to be variable following a normal distribution, so as to better approach the real situation where there is no guarantee of maintaining a perfect contact during the mapping. Our results show that reducing misalignment among u-EGMs within the clique improves fibrosis detection ability of the proposed EIGDR-based index. This is in agreement with other studies where time alignment of b-EGMs has shown to be beneficial for electroanatomical mapping strategies robustness [8]. The index R A provides comparable fibrosis detection accuracy to the one of maximum bipolar voltage maps when u-EGMs are not affected by noise, and better when high noise levels are present ( v ≥ 23.2 V ), for both fixed and variable electrode-to-tissue distances.
Results obtained by considering the three MEA orientations separately reinforce the consideration of R A as an index worth to be analyzed in extended studies with real recordings for discriminating between fibrotic and normal areas, also pointing out the larger impact of catheter orientation in bipolar amplitudes than in EIGDR metrics.
In addition, if the evaluation is performed without exclusion of cliques which have electrodes inside and outside the fibrotic patch, similar conclusions, with reduced difference ranges, can be drawn. The ratio R A still shows higher ACC in the largest noise contamination and with variable electrodeto-tissue distance. Moreover, EIGDR-based markers reveal to be more robust to the effect of variable distance than bipolar maps, especially when u-EGMs are not affected by noise.
Regarding threshold values corresponding to the maximum fibrosis detection accuracy, our findings reveal that the thresholds needed to maximize accuracy of bipolar maps V b-y and V b-m are greater than the one typically used in clinical settings (0.5 V ), whereas V b-x presents lower voltage threshold. This is explained by the fact that there is no projection of wavefront propagation along the x-axis of the MEA when propagation orientation is Ψ = 0 • .

Performance evaluation of fibrosis markers with clinical data
In the present study, we also tested the ability of the EIGDRbased markers to characterize the fibrotic substrate considering different mapping points acquired with a PentaRay ® catheter within fibrotic and non-fibrotic areas over the LA anatomical map, using MRI as reference for that purpose. The electrode clique organization, referred to a fixed structure catheter like the Advisor TM HD Grid, was extended to a flexible structure catheter like the PentaRay ® , where the inter-electrode spacing within the clique may vary at different acquisition points. Nevertheless, this does not represent a problem for the proposed method, as it is not dependent on electrode orientation.
Preliminary findings obtained from real u-EGMs in this paper reveal that the ratio based on time-aligned u-EGMs R A is the only EIGDR-based marker between F and NF mapping points, also showing better discrimination power than bipolar amplitudes V b and V b-m typically used in clinical settings. In addition, R A , together with ΔR A , proved to be the only indices capable to globally discriminate fibrosis from non-fibrotic tissue, regardless of the mapping points and cliques/bipoles considered at each of them.

Limitations
Several limitations of this study need to be highlighted. First, we simulated a single scenario reproducing a simple propagation pattern in a 2D atrial model, which largely simplifies the real 3D anatomical and electrophysiological situations.
Although a single plane wavefront that propagates in a homogeneous tissue lends itself well to approximating the propagation during pacing, the previous considerations do not allow us to extend quantitatively the results to other propagation patterns, such as circular waves, wave collisions, reentrant wave fronts, among others, and model conditions, including conduction anisotropy or patchy fibrosis. Nevertheless, even #2 #5 Fig. 6 3D reconstruction of the LA geometry (gray mesh) and corresponding co-registered MRI, showing the different regional distribution patterns of gadolinium (red areas: latest contrast enhancement, blue areas: absence of latest contrast enhancement). In the geometrical mesh, two of the mapping points acquired and considered in the analysis (point #5 at fibrosis, point #2 at non-fibrosis) are marked and color-coded according to their corresponding bipolar peak-to-peak amplitude. For each of them, the atrial activation windows extracted from the twenty filtered u-EGMs recorded with the PentaRay ® catheter are also displayed. Note that not all displayed u-EGMs recorded at a particular catheter site belong to a clique, see Section 2.3, and therefore affect the EIGDR indices and bipolar amplitude computations 13 × 13 ground-truth masks for evaluating fibrosis detection ability of maps performed with 2 × 2 and 3 × 3 cliques, respectively. Green squares represent the pixels corresponding to cliques with some electrodes inside and some outside the fibrotic patch, i.e. those cliques lying in the border separating the fibrotic patch from non-fibrotic tissue, which were excluded from the evaluation if it is well-established that the underlying propagation pattern strongly influences EGMs morphology and their spatiotemporal information, we expect that it does not largely affect the local EIGDR computation. This is because we hypothesize that the correlation between the presence of fibrosis and the morphology dispersion of signals in electrode cliques is well modeled by a waveform assumed to be locally plane and homogeneous, irrespective of the global waveform distribution across the complete tissue. For the same reasons, we considered that the EIGDR approach would not be largely affected by the shape and size of the fibrotic patch. Note that the proposed intraclique time alignment of the u-EGMs compensates the effect of different u-EGM arrival times on EIGDR. This leaves EIGDR to mostly represent spatial relationships differences among u-EGMs within each clique. In this work, only the effect of broad-band noise affecting u-EGMs was considered, while specific periodic types of noise were not considered. It must be noted that far-field disturbances due to ventricular depolarization did not occur during atrial activation in sinus rhythm.
Lastly, results presented with real signals represent a proof of concept, but increased sample size need to be considered in order to elucidate whether the use of the EIGDR-based approach is advantageous for fibrosis detection in clinical settings.

Conclusions
In this paper we demonstrated that mapping strategies based on the EIGDR method are able to discriminate fibrotic from non-fibrotic tissue. In simulation, they attain comparable performance to map obtained by combining the b-EGMs amplitudes along the two directions of the MEA for low noise levels, when assuming both fixed and variable distance between the electrode grid and the tissue. Nevertheless, they outperform bipolar maps when higher noise levels are added. Moreover, performance of 3 × 3 electrode cliques outperforms the 2 × 2 cliques one and fibrosis detection benefits from the previous time alignment of u-EGMs. With clinical data, EIGDR approach showed promising results in discriminating fibrotic and non-fibrotic mapping points, especially when u-EGMs are previously aligned in time. Both scenarios studied lead to choose R A as EIGDR biomarker for fibrosis discrimination.