Unipolar Electrogram Eigenvalue Distribution Analysis for the Identification of Atrial Fibrosis

Unipolar Electrogram Eigenvalue Distribution Analysis for the Identification of Atrial Fibrosis Jennifer Riccio1, Sara Rocher2, Laura Mart´ ınez-Mateu2, Alejandro Alcaine3,1, Javier Saiz2, Juan Pablo Mart´ ınez1,3, Pablo Laguna1,3 1BSICoS, I3A, IIS Arag´ on, Universidad de Zaragoza, Zaragoza, Spain 2Ci2B, Universitat Polit` ecnica de Val` encia, Valencia, Spain 3CIBER en Bioingenier´ ıa, Biomateriales y Nanomedicina (CIBER-BBN), Spain Abstract Atrial fibrosis plays an important role in the pathogenesis of atrial fibrillation (AF). Low bipolar electrograms (bEGMs) peak-to-peak voltage areas indicate scar tissue and are considered targets for AF substrate ablation. However, this approach ignores the spatiotemporal information embedded in the signal and the dependence of b-EGMs on catheter orientation. This work proposes an approach to detect fibrosis based on the eigenvalue dominance ratio (EIGDR) in an ensemble (clique) of unipolar electrograms (u-EGMs). A 2-D tissue with a central circular patch of fibrosis has been simulated using the Courtemanche cellular model. Maps of EIGDR have been computed using two sizes of electrode cliques, from the original u-EGMs within the ensemble or after a time alignment of these signals. Performance of each map in detecting fibrosis has been evaluated using receiver operating characteristic curves and detection accuracy. Best results achieve an area under the curve (AUC) of 0.98 and an accuracy (ACC) of 1 when we use as marker the gain in eigenvalue dominance produced by the ensemble alignment. 1. Introduction Atrial fibrosis represents a structural abnormality of the atrium, which alters the electrical conduction and excitability of the tissue. Fibroblasts proliferation and their secretion of extracellular matrix proteins, such as collagen, characterize fibrotic tissue [1]. These proteins are mainly involved in the reparative process to replace damaged myocardial parenchyma [2]. Atrial fibrosis is observed to be closely related to atrial fibrillation (AF), even if their causal relationship is still challenging [3]. On the one hand, an extensive fibrotic process in the atrium can promote persistent AF [3]; on the other hand, structural atrial remodeling found in AF produces fibrosis and alters tissue function [2]. Electrophysiologically, atrial fibrosis produces low-voltage intracardiac electrograms (EGMs), which can be identified using electroanatomical mapping (EAM) [4]. Peak-to-peak bipolar voltage maps can be constructed through data obtained during substrate mapping [5], being bipolar voltage an interesting marker during sinus rhythm (SR) as well as in AF. Low-voltage areas are typically defined as those with peak-to-peak bipolar voltage below 0.5 mV during SR. However, some drawbacks to this procedure should be put in evidence. First, peak-to-peak voltage measure does not provide information about morphological features or temporal trend embedded in the signal; therefore, voltage thresholding does not take into account the presence of underlying abnormalities in the atria. Second, the methodology to define low voltage areas has not been standardised [6]. Third, spatial heterogeneity is not accounted for. Fourth, low bipolar voltage can also be influenced by other factors than fibrosis, such as activation direction, electrode size, interelectrode distance and filtering [4], [6], as well as by technical problems such as poor electrode contact in anatomically difficult sites (e.g. the pulmonary veins) or its instability in time. In order to overcome these limitations, in this work we propose the dominant-to-remaining eigenvalue dominance ratio (EIGDR) of unipolar EGMs (u-EGMs) as a measure of the voltage wavefront roughness and correlate it with the presence of fibrosis. We compute maps of eigenvalues ratios considering two different electrode clique arrangements (ECA) and we evaluated the ability of each map to detect a fibrosis patch in the context of a simulation study. Maps have been created from the whole length of u-EGMs, using 0oand 45ocatheter orientations with respect to the wavefront propagation. 2. Materials A 2-D atrial tissue of 4x4 cm of hexahedric elements simulated with 100 µm resolution using the Courtemanche Computing in Cardiology 2020; Vol 47 Page 1 ISSN: 2325-887X DOI: 10.22489/CinC.2020.434 cellular model [7] has been used. Conduction heterogeneity induced by chronic AF in left atrium has been considered in the model. The simulated tissue includes a circular patch of diffuse fibrosis having a diameter of 2 cm, where 20% nodes have randomly been assigned the Maleckar model for fibroblasts [8] and a conductivity reduction of 30%. Let xk(n)be the u-EGMs computed with a sampling frequency of 1 kHz at a high-density multi-electrode array (MEA) of 15×15 electrodes, k∈ {1,...,225}, located at sites (i, j), i, j ∈ {1, .., 15}. The MEA has an interelectrode distance d= 2 mm, is centered in the tissue slice and located at 1 mm distance from the tissue surface. Each simulated u-EGM is 500 ms long and contains a single activation (depolarization plus repolarization) corresponding to one sinus beat. 3. Methods 3.1. Eigenvalue analysis In this work, we propose and assess EIGDR from uEGMs at each 3×3 and 2×2 clique from the MEA, as fibrosis markers. Eigenvalues were obtained from the spatial covariance matrix of both original and aligned uEGMs within clique ensembles. Signals have been aligned as proposed in [9], according to the maximum crosscorrelation with respect to the highest amplitude u-EGM. For each electrodes clique, the ratio Rof the dominant-toremaining eigenvalues, computed as: R=λ1 PK k=2 λk ,(1) where Kis the number of u-EGMs xk(n)in the clique ensemble (4 or 9 in this study), is estimated to quantify EIGDR. For a theoretical analysis of EIGDR, the model xk(n) = αks(n−τk) + fk(n) + vk(n)(2) is considered for the u-EGMs xk(n)at the kth electrode, where s(n)is the clean and space invariant u-EGM in the case of a plane wave propagation, τkis the delay of the kth u-EGM s(n−τk)with respect to the time reference in the ensemble, αkis a parameter accounting for u-EGM amplitude (decreased at fibrosis, αk<1, with respect to normal tissue, αk= 1), fk(n)represents the u-EGM fibrotic component (absent in normal tissue), and vk(n)is the noise component of the kth u-EGM. The energy of s(n)is denoted by Es, and Es0represents the energy of its derivative s0(n). The delays τkare characterized by their variance, β2σ2 θ, with β > 1 in fibrosis, a factor accounting inversely for fibrosis generated speed reduction relative to normal tissue where β=1. αkis modelled as random variable with mean E[αk] = αand variance σ2 α. Noise is considered to be zero-mean, Gaussian, white and uncorrelated with τkand fk(n), with variance σ2 v. The variance of the zero-mean fibrotic component across clique electrodes is denoted by σ2 f. Four u-EGMs scenarios for non-aligned (NA) and aligned (A) u-EGMs at non-fibrotic (NF) and fibrotic (F) areas are considered and their approximate theoretical eigenvalues are derived according to the procedure presented in [10]. Table 1 shows the eigenvalues λkand EIGDR for the four scenarios. Fibrosis markers: •From the derivations in equations presented in Table 1, it is clear that R>RF, as result of three concomitant effects appearing simultaneously at fibrosis: a) higher morphology dispersion, σ2 f, b) lower amplitude, α < 1, and c) larger delays resulting in larger misalignment dispersion β > 1. Then Ris proposed as one marker for fibrosis detection. •Similarly, for the same two first reasons, it also results that RA>RA F, suggesting RAas other fibrosis detection marker. •Analyzing the ratio ∆RFbetween EIGDR in nonfibrotic areas with respect to fibrotic ones, representing the eigenvalue concentration lost by fibrosis, and then a measure of the separability power of RFas a marker for fibrosis, we obtain for misaligned u-EGMs, σ2 θ>0: ∆RF=R RF ≈β2σ2 θEs0+N(σ2 v+σ2 f) (α2+σ2 α) σ2 θEs0+Nσ2 v .(3) To evaluate how ∆RFvaries with the level of misalignment σ2 θwe perform the derivative, obtaining: ∂∆RF ∂σ2 θ ≈ −Es0Nσ2 f+σ2 v1−β2α2+σ2 α α2+σ2 α(σ2 θEs0+Nσ2 v)2. (4) The term 1−β2α2+σ2 α is typically >0since in fibrosis βcan get values up to 2 and αup to 1/8 [11], resulting that ∂∆RF ∂σ2 θ <0, justifying the advantage of alignment, since the higher the misalignment σ2 θ, the lower the separability capacity of RFto discriminate between fibrosis and non-fibrosis, suggesting that RA Fis better suited marker than RF. •Alternatively, the ratio ∆RAbetween EIGDR before and after alignment, representing the gain in eigenvalue concentration produced by the ensemble alignment, is considered. In the case of fibrosis, it is reached the value: ∆RA=RA F RF ≈ Esβ2σ2 θα2+σ2 αEs0+N(σ2 v+σ2 f) N(σ2 v+σ2 f)(Es−β2σ2 θEs0) . (5) To evaluate how ∆RAvaries with the level of fibrosis σ2 f we also perform the derivative, obtaining: Page 2 u-EGM model λkEIGDR NA, NF xk(n) = s(n−τk) + vk(n)λk≈   (Es−σ2 θEs0)K/N +σ2 v, k = 1; σ2 θEs0K/N +σ2 v, k = 2; σ2 v, k = 3,...,K, R ≈ Es−σ2 θEs0 σ2 θEs0+Nσ2 v . A, NF xk(n) = s(n) + vk(n)λk≈EsK/N +σ2 v, k = 1; σ2 v, k = 2,...,K, RA≈Es Nσ2 v. NA, F xk(n)=αks(n−τk)+fk(n)+vk(n)λk≈     α2+σ2 α(Es−β2σ2 θEs0)K/N +σ2 v+σ2 f, k = 1; α2+σ2 αβ2σ2 θEs0K/N +σ2 v+σ2 f, k = 2; σ2 v+σ2 f, k = 3, ., K RF≈Es−β2σ2 θEs0 β2σ2 θEs0+ N(σ2 v+σ2 f) (α2+σ2 α) A, F xk(n) = αks(n) + fk(n) + vk(n)λk≈(α2+σ2 αEsK/N +σ2 v+σ2 f, k = 1; σ2 v+σ2 f, k = 2,...,K, RA F≈Es N(σ2 v+σ2 f) (α2+σ2 α) . Table 1: Models for non-aligned (NA) and aligned (A) u-EGMs at non-fibrotic (NF) and fibrotic (F) areas, with their respective eigenvalues λkand eigenvalue dominance ratios EIGDR computed following the procedure presented in [10] ∂∆RA ∂σ2 f ≈−EsNβ2σ2 θα2+σ2 αEs0 N(σ2 v+σ2 f)2(Es−β2σ2 θEs0) .(6) Since for small delays τk,Esβ2σ2 θEs0it results that ∂∆RA ∂σ2 f <0, making this marker ∆RAbecoming smaller the higher the fibrotic component σ2 f, justifying to consider it as a potential fibrosis marker. Also, ∂∆RA ∂α2>0, so the larger the fibrosis, (reduced α2), the further gets ∆RAreduced. However, ∂∆RA ∂β2>0, and since βincreases in fibrosis, it results in a counteracting effect for the marker sensitivity to fibrosis. Since fibrosis effects on u-EGM amplitude, α2, and morphology, σ2 f, are much more marked than on conduction velocity, β2, [11], it is expected that the first two tendencies dominate, making the marker ∆RAto reduce with fibrosis. Simulation and real experiments can elucidate the behaviour in practice. 3.2. Assessment of R,RAand ∆RAfor fibrosis detection Maps of R,RAand ∆RAhave been created processing the complete MEA with the electrode lines in two orientations with respect to the waveform propagation direction, parallel (0o) and oblique (45o) and two ECA (3×3 and 2×2). The 3×3ECA provides one EIGDR for each squared group of nine electrodes with diagonal vertices at (i, j), and (i+2, j+2), i, j ∈ {1,...,13}, giving a total of 13 ×13 pixel maps of markers R,RAor ∆RA. The 2×2 ECA provides values at each squared group of four electrodes with diagonal vertices at (i, j)and (i+ 1, j + 1), resulting in maps of 14 ×14 pixels, i, j ∈ {1,...,14}. Receiver operating characteristic (ROC) curves have been used to evaluate the markers ability in discriminating fibrotic from non-fibrotic areas. For that purpose, a ground truth mask has been created labelling these areas. Cliques lying in the tissue interface, i.e., in the line separating the fibrotic patch from non-fibrotic tissue, were not considered in the evaluation. For each map, thresholds for fibrosis identification have been varied to compute the ROC curve. The area under the curve (AUC) and the maximum accuracy (ACC) have been computed for each map, as a measure of the ability of the marker to detect fibrosis. AUC/ACC Clique Angle R RA∆RA 3×30o0.88/0.86 0.96/0.95 0.97/0.95 45o0.77/0.75 0.96/0.93 0.98/1 2×20o0.84/0.80 0.85/0.79 0.72/0.70 45o0.76/0.71 0.96/0.93 0.94/0.89 Table 2: AUC and ACC of the three different metrics 4. Results Figure 1 shows the mean and standard deviation (SD) for each studied marker R,RAand ∆RAat fibrotic and healthy simulated tissues, using 2×2 or 3×3 cliques. The three markers were significantly lower in the healthy than in the fibrotic tissue (Wilcoxon rank-sum test, p < 0.05). RApresents the higher ratio of fibrotic to healthy tissue γ= 2.66 (2.51) for 3×3 (2×2) cliques, being ∆RA the marker with more significant differences between both types of tissue (p < 0.05). Table 2 contains AUC and ACC for the markers considered. Results show higher values when time alignment of u-EGMs is performed, especially when MEA has an orientation of 45owith respect to the wavefront direction. The best discrimination power is obtained by ∆RA, with AUC = 0.97 (0.98) for parallel (diagonal, where ACC = 1) orientation, using 3×3ECA. Figure 2 shows the maps of R, RAand ∆RAobtained with 3×3 ECA and with parallel catheter orientation. In the lower panels, the detected fibrotic areas are shown, using the thresholds that maximize the detection accuracy of each marker. 5. Discussion and Conclusions In this work, the eigenvalue dominance ratios, EIGDR, of u-EGMs have been proposed and evaluated using a simPage 3 3x3 2x2 0 5 10 15 20 F NF 3x3 2x2 0 100 200 300 400 F NF 3x3 2x2 0 5 10 15 20 25 F NF = 1.72 = 1.56 = 2.51 = 2.66 = 1.49 = 1.55 R RA∆RA Figure 1: EIGDR markers (R,RAand ∆RA) mean and SD from the two cliques, 3×3 and 2×2, considering both MEA orientations, 0oand 45o.γis the factor relating the mean from fibrotic to non-fibrotic areas for each marker and clique size. 0 1 2 3 4 5 6 7 8 9 10111213 0 1 2 3 4 5 6 7 8 9 10 11 12 13 4 5 6 7 8 0 1 2 3 4 5 6 7 8 9 10111213 0 1 2 3 4 5 6 7 8 9 10 11 12 13 40 60 80 100 120 140 0 1 2 3 4 5 6 7 8 9 10111213 0 1 2 3 4 5 6 7 8 9 10 11 12 13 8 10 12 14 16 18 R RA∆RA 0 1 2 3 4 5 6 7 8 9 10111213 0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 10111213 0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7 8 9 10111213 0 1 2 3 4 5 6 7 8 9 10 11 12 13 0 0.2 0.4 0.6 0.8 1 Figure 2: Top panels: maps of R,RAand ∆RAfrom 3×3 cliques for catheter orientation of 0o. Circle encompasses fibrotic area. Lower panels: detected fibrotic areas (brown), using the thresholds that maximize detection accuracy of each marker. Blue (brown) color inside the circle denotes FN (TP), while outside denotes TN (FP) detection, respectively. ulated atrial tissue including areas with fibrosis and others of normal tissue, in order to discriminate them. In clinical setting, bipolar voltage is commonly used as a surrogate of atrial fibrosis, but the phenomenon is much more complex and voltage cannot be considered as a substitute of fibrosis as assessed by MRI. Here, three different EIGDR markers have been studied, evaluating their ability in two possible ECA for two catheter-to-wavefront orientations and considering depolarization and repolarization of the u-EGMs. The proposed markers are good candidates for detecting fibrotic tissue. Results in terms of AUC and ACC confirm the hypothesis that reducing misalignment is beneficial and suggest ∆RA, representing the eigenvalue dominance gain by alignment, as the better suited for fibrotic areas identification, especially when the 3×3ECA is used. However, these results need to be complemented with other simulation configurations, such as smaller fibrotic areas and patchy fibrosis, as well as with real data, where different values for the electrodes size can be tested and the quality of tissue-electrode contact could also be considered. Acknowledgments Funding comes from EU Programme H2020 under the Marie Sklodowska-Curie Grant No 766082 (MY-ATRIA), Gobierno de Arag´ on (BSICoS Group T39-20R) cofunded by FEDER 2014-2020 “Building Europe from Aragon”, fellowship ACIF/2018/174 from Generalitat Valenciana, and PID2019-104881RB-I00 from MICINN, Spain. References [1] Tzeis S, Asvestas D, Vardas P. Atrial fibrosis: translational considerations for the management of AF patients. AER Journal 2019; 8(1):37–41. [2] Burstein B, Nattel S. Atrial fibrosis: mechanisms and clinical relevance in atrial fibrillation. JACC 2008; 51(8):802–809. [3] Platonov PG. Atrial fibrosis: an obligatory component of arrhythmia mechanisms in atrial fibrillation?. J Geriatr Cardiol. 2017; 14(4):233–237. [4] Rodr´ ıguez-Ma˜ nero M et al. Validating left atrial low voltage areas during atrial fibrillation and atrial flutter using multielectrode automated electroanatomic mapping. JACC: Clinical Electrophysiology 2018; 4(12):1541-1552. [5] Bhakta D, Miller JM. Principles of electroanatomic mapping. Indian Pacing Electrophysiol J. 2008; 8(1):32-50. [6] Sim I, Bishop M, O’Neill M, Williams SE. Left atrial voltage mapping: defining and targeting the atrial fibrillation substrate. J Interv Card Electrophysiol 2019; 56:213-227. [7] Courtemanche M, Ramirez RJ, Nattel S. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. Am. J. Physiol. 1998; 275(1), H301-21. [8] Maleckar MM, Greenstein JL, Giles WR, Trayanova NA. Electrotonic coupling between human atrial myocytes and fibroblasts alters myocyte excitability and repolarization. Biophys. J. 2009; 97:2179–2190. [9] S¨ ornmo L, Laguna P. Bioelectrical Signal Processing in Cardiac and Neurological Applications, Elsevier (Academic Press), Amsterdam, 2005. [10] Laguna P et al. Eigenvalue-based time delay estimation of repetitive biomedical signals. Digital Signal Processing 2018; 75:107-119. [11] Vigmond E et al. Percolation as a mechanism to explain atrial fractionated electrograms and reentry in a fibrosis model based on imaging data. Heart Rhythm 2016; 13(7):1536-1543. Address for correspondence: Jennifer Riccio C/ Mariano Esquillor s/n, Edificio I+D+i, Lab 6.1.01 50018 Zaragoza, Spain. [email protected] Page 4