Characterization of Atrial Propagation Patterns and Fibrotic Substrate With a Modified Omnipolar Electrogram Strategy in Multi-Electrode Arrays

Introduction: The omnipolar electrogram method was recently proposed to try to generate orientation-independent electrograms. It estimates the electric field from the bipolar electrograms of a clique, under the assumption of locally plane and homogeneous propagation. The local electric field evolution over time describes a loop trajectory from which omnipolar signals in the propagation direction, substrate and propagation features, are derived. In this work, we propose substrate and conduction velocity mapping modalities based on a modified version of the omnipolar electrogram method, which aims to reduce orientation-dependent residual components in the standard approach. Methods: A simulated electrical propagation in 2D, with a tissue including a circular patch of diffuse fibrosis, was used for validation. Unipolar electrograms were calculated in a multi-electrode array, also deriving bipolar electrograms along the two main directions of the grid. Simulated bipolar electrograms were also contaminated with real noise, to assess the robustness of the mapping strategies against noise. The performance of the maps in identifying fibrosis and in reproducing unipolar reference voltage maps was evaluated. Bipolar voltage maps were also considered for performance comparison. Results: Results show that the modified omnipolar mapping strategies are more accurate and robust against noise than bipolar and standard omnipolar maps in fibrosis detection (accuracies higher than 85 vs. 80% and 70%, respectively). They present better correlation with unipolar reference voltage maps than bipolar and original omnipolar maps (Pearson's correlations higher than 0.75 vs. 0.60 and 0.70, respectively). Conclusion: The modified omnipolar method improves fibrosis detection, characterization of substrate and propagation, also reducing the residual sensitivity to directionality over the standard approach and improving robustness against noise. Nevertheless, studies with real electrograms will elucidate its impact in catheter ablation interventions.


INTRODUCTION
Diagnosis and treatment of a wide range of atrial arrhythmias, including atrial fibrillation, require extraction of substrate and propagation features from intracardiac electrograms (EGMs) (Houben and Alessie, 2006). Mechanisms underlying atrial fibrillation are not yet completely clear (Calkins et al., 2017) and no single-one can comprehensively explain the arrhythmia. The electrical isolation of focal trigger sites in the pulmonary veins (Keane and Ruskin, 2002) may abrogate early stage atrial fibrillation when anti-arrhythmic drug therapy becomes ineffective (Wazni et al., 2005). However, in more advanced stages of persistent atrial fibrillation, additional arrythmogenic substrate and mechanisms such as fibrosis are likely present (Lau et al., 2017), ablation of which was shown to improve success rates in patients (Jadidi et al., 2016).
Due to its unstable nature from beat to beat, atrial fibrillation is better investigated through simultaneous EGMs, recorded by high-density multi-electrode catheters from multiple sites at the same time, rather than by analyzing sequential EGMs, each recorded at a different time within the mapping procedure. From such signals, meaningful features for characterizing atrial substrate and propagation pattern are extracted and mapped through 3D electroanatomic mapping systems, thus helping physicians to visualize the locations that are generating the erratic propagation.
Characterization of atrial substrate is typically performed by mapping the peak-to-peak amplitude of bipolar EGMs (b-EGMs). This strategy identifies as fibrosis  those areas having voltage lower than 0.5 mV (Sim et al., 2019). The fibrotic areas determined in this way constitute a potential substrate for atrial fibrillation maintenance (Burstein and Nattel, 2008) and, consequently, a target for ablation, in combination with pulmonary veins isolation (Haldar et al., 2017). However, the dependence of b-EGMs on multiple factors, including catheter orientation with respect to wavefront direction, electrode size, inter-electrode spacing and contact with the tissue, affect bipolar amplitude values . Moreover, the proposed threshold to define low-voltage areas has not received histological or electrophysiological evidence (Yamaguchi et al., 2019).
In current practice, characterization of atrial propagation patterns is performed by assessing local activation times (Magtibay et al., 2019), using unipolar EGMs (u-EGMs) or b-EGMs. However, both u-EGMs and b-EGMs present limitations because: (1) u-EGMs are sensitive to electric far-field disturbances due to other large cardiac structures such as the ventricles; (2) b-EGMs depend on the angle between the catheter and wavefront propagation direction; (3) both are sensitive to local recording noise.
The recently proposed Omnipolar EGM (OP-EGM) method tries to cope with the problems affecting unipolar and bipolar measurements (Deno et al., 2017). It is based on the estimation of the electric field from all the b-EGMs locally recorded at each group of nearby electrodes, referred as clique by Deno et al. (2017), under the assumption of locally plane and homogeneous propagation within the clique. Two types of cliques are used in this work: (a) square cliques, which consist of four-electrode sets forming a square, and (b) triangular cliques, which consist of three-electrode sets formed by an electrode and two adjacent ones, one in each direction of the 2D plane. From these electrode sets, omnipolar signals and parameters are derived and associated to a single point virtually located at the center of the clique. The electric field evolution over time within each clique describes a loop which depends on the propagation direction. OP-EGM signals are obtained as projections of this loop onto the propagation direction, as if a virtual bipole in that direction had been used. In this way, the dependence on the catheter orientation is theoretically reduced. Substrate and propagation parameters such as peak-to-peak amplitude of OP-EGM signals and conduction velocity (CV) are also estimated at the clique level without requiring previous estimation of local activation times.
In this work, we proposed voltage and CV mapping strategies based on modifications of the OP-EGM method (MOP-EGM) by using two clique configurations: square and triangular cliques. These modifications are hypothesized to improve the accuracy and robustness of the voltage and CV estimates and to reduce the error induced by the b-EGMs dependence on catheter orientation, which is not fully compensated by the OP-EGM approach. Specifically, the aims of this study are: (1) to assess the estimates of CV and propagation direction within each clique with the proposed MOP-EGM method as compare to the ones obtained by the original OP-EGM; (2) to propose mapping strategies based on omnipolar EGMs estimations to characterize the atrial substrate; and (3) to assess the ability of the proposed maps in detecting a fibrotic patch and in reproducing unipolar voltage maps.

2D Atrial Sheet Model
In order to assess the proposed mapping strategies, a 2D atrial tissue of 4×4 cm under chronic atrial fibrillation (cAF) conditions was simulated 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. Inside this circular patch, a diffuse fibrosis pattern was randomly defined following a uniform distribution. All the non-fibrotic nodes were assigned the Courtemanche cellular model stabilized during 1 min at basic cycle length (BCL) 500 ms (Courtemanche et al., 1998). To reproduce left atrial action potential (AP) morphology, the maximum ionic conductance of the rapid delayed rectifier potassium current (I Kr ) was selected 1.6-fold greater than the original right atrial model (Li et al., 2001;Martinez-Mateu et al., 2018). Additionally, electrical remodeling induced by cAF 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 Supplementary Table 1), as in previous computational studies (Tobón et al., 2013;Martinez-Mateu et al., 2018). Twenty percent of the nodes within the fibrotic patch were assigned the Maleckar fibroblast model (Maleckar et al., 2009). The intercellular conduction velocity was decreased to 30% in all elements with at least one fibroblast node . Although the percentage of atrial fibrosis is very patient-dependent, the 20% represents the threshold value between stage II and stage III according to the Utah classification (Chelu et al., 2018), therefore being a realistic percentage for this study. Supplementary Figure 1 shows the effect of adding the cAF electrical remodeling and the fibroblast coupling in the Courtemanche model. The APs were registered in two different myocytes from the mesh, one outside the fibrotic patch and the other inside the fibrotic patch coupled with two fibroblasts. Electrical remodeling produces a 55% reduction in duration measured at 90% repolarization (APD 90 ) (248 vs. 111 ms), in concordance with experimental data (Bosch et al., 1999). Fibroblast coupling with cAF myocytes makes resting potential less negative (83 vs. 78 mV) and elongates the APD 90 (111 vs. 120 ms). Simulations were run in ELVIRA software (Heidenreich et al., 2010). The monodomain formulation was solved using the operator splitting numerical scheme with a constant time step dt = 0.01 ms and a spatial resolution dx = 0.1 mm.
Recording electrodes were distributed over the simulated tissue, mimicking a N × N high-density multi-electrode array (MEA) (N = 15), where the inter-electrode distance was d = 2 mm. The simulated tissue geometry, including the MEA location and the fibrotic patch, is shown in Figure 1A. Figure 1B depicts a generic square clique. To keep the same simplified notation in each position within the MEA, its lower left electrode, corresponding to location (i, j) within the MEA, is referred to as electrode 1. The rest electrodes of the clique are numbered from left to right and bottom to top. Therefore, electrodes 2, 3, and 4 correspond to locations (i + 1, j), (i, j + 1), and (i + 1, j + 1), respectively.
The simulated electrode grid is located so that its central electrode corresponds to the center of both the tissue slice and the fibrotic patch, and the MEA is rotated by an angle with respect to the tissue fiber orientation. Three different MEA orientations with respect to the tissue preferential direction were considered: = 0 • , = 30 • , and = 45 • , which are illustrated in Figure 2A. For the case when = 0 • , the propagating wavefront direction was depicted in Figure 2B, leftmost, where the fibrotic patch at the center of the tissue is also marked with a diffuse colored area. Moreover, for the same case, the propagation pattern was shown in Figure 2B, rightmost, by means of the map of local activation times, estimated as the maximal negative slope of the available synthetic u-EGMs (Paul et al., 1990). A video content showing propagating wavefront can also be found in the Supplementary Material.
The MEA grid was placed at z = 1 mm distance from the tissue surface. This represents a plausible compromise between recording realistic amplitude EGMs and the real clinical setting situation where there is no guarantee of maintaining a perfect electrode-tissue contact during the mapping.

Synthetic Signals
Unipolar EGMs u i,j (t) were simulated as generated by the propagation wave when passing by electrodes located at sites (i, j), (i ∈ {1, ..., N}, j ∈ {1, ..., N}) of the grid. They were computed in a volumetric tissue-blood model with a temporal resolution of 1 ms, as in Martinez-Mateu et al. (2019), by using an approximation of the bidomain formulation in two steps (Keller et al., 2010;Martinez-Mateu et al., 2018). By assuming equal anisotropy ratios for the intracellular and extracellular conductance tensors, bidomain equations can be decoupled accounting one of them for changes in the transmembrane potential and the other one for changes in the extracellular potential. Therefore, in a first step, transmembrane potential was solved by the monodomain approach. Then, in a second step, the already calculated transmembrane potential was used to obtain the extracellular potential, considering the tissue immersed in a non-conductive bath.
Simulated signals were obtained with a sampling frequency of 1 kHz, and had a duration of 0.5 s, containing one single activation (depolarization plus repolarization). Bipolar EGMs b x i,j (t) and b y i,j (t) were derived along the two MEA directions as follows: Bipolar signals, as defined in Equation (1), are used to compute bipolar voltage maps for performance comparison, as described in subsections 2.12, 2.13.
The omnipolar methodological developments assume a locally plane wave propagation within a clique, with a velocity vector v = v − → u v , where − → u v denotes a unit direction vector whose spatial dependence has been omitted for notation simplicity. In the 2D Cartesian coordinate system defined by the unit vectors in the main directions of the MEA ( − → u x and − → u y ), the local plane wave propagation direction within the clique having electrode (i, j) in its lower left corner,  − → u v , forms an angle θ (i,j) with the direction of − → u y , as illustrated in Figure 2C. It must be noted that the angle is unique for each 2D catheter configuration and should not be mistaken for θ (i,j) , which depends on local propagation at (i, j) electrode position, see black arrows direction variations in Figure 2B, leftmost.
Frontiers in Physiology | www.frontiersin.org Examples of unipolar and bipolar EGMs used in this study, corresponding to the case when the global propagation direction of the wave and MEA direction form an angle = 0 • , are shown in Supplementary Figure 2. In that case, b-EGMs are reported with added noise as well (see also subsection 2.11).

Activation Times and Relation Between Electrograms Under Plane Wave Assumption
Let us consider any four-electrode clique within the MEA, like the one in Figure 1B. Let t c be the reference time at which the center c of the selected square is activated. Assuming constant velocity within the arrangement, the activation times at the four electrodes are: Let φ(t) be the unipolar voltage waveform generated by a plane wave when passing by a given position at time t = 0. When one or more electrodes within the square of Figure 1B are activated by φ(t) passage, the following u-EGMs and b-EGMs are modeled as: It should be noted that bipolar signals in Equation (4) are derived within each clique along x direction (b 12 (t) and b 34 (t)), y direction (b 13 (t) and b 24 (t)), as well as along the directions of its diagonals (b 14 (t) and b 23 (t)). Depending on the waveform arrival time to each electrode t m , the activation time of each b-EGM b mn (t) corresponds to the passing of the wave by the middle point between electrodes m and n. Therefore, relative delays (misalignment) occur between the different b mn (t). By using the Taylor's series expansion of u m (t) = φ(t − t m ) around t m = t c , we can approximate: where t mn = t n −t m are the differences between activation times at electrodes m and n, m, n ∈ {1, 2, 3, 4}, and t m = t m − t c are the differences between activation time at electrode m and the reference activation time of the center of the square. By replacing Equation (2) in Equation (6) and operating, the following b-EGMs expressions are derived: From these expressions we can relate the second order approximations of the different pairs of b-EGMs b mn (t), showing the previously mentioned delay between the different bipolar signals. If we relate b 12 (t) and b 13 (t), measured in different x and y directions, we obtain: Equation (13) shows that the b-EGM along y-direction (1-3) is a scaled version of the b-EGM along x-direction (1-2) by a factor cos θ sin θ and delayed by a time . This corresponds to the delay between the pass of the wave by the middle points of each electrode pair, separated d/ √ 2. As expected, for θ = 90 • (propagation in x-direction), b 13 (t) = 0, while for θ = 45 • , the delay is τ = 0, meaning that both bipolar EGMs in the x-and y-directions are activated at the same time. Similarly, if we relate b 12 (t) with b 34 (t), both along the x-direction, we obtain: Equation (14) shows that b-EGMs along parallel directions have the same amplitudes (as expected), and the delay τ = d v cos θ Frontiers in Physiology | www.frontiersin.org ranges from 0, when the propagation is in the same direction than the electrode pairs (θ = 90 • ), which are therefore activated at the same time, to a maximum of d/v, when the propagation is orthogonal to the direction of the electrode pairs (θ = 0 • ).

OP-EGM Framework
The purpose of the OP-EGM method is to obtain EGM amplitude and propagation features which are invariant to the relative angle between the propagation direction and the catheter (Deno et al., 2017). This method is based on the relationship between the spatial gradient of voltage φ(t), and the electric field (E-field) at the extracellular-myocardial interface. Let us consider a locally plane myocardial surface, defined by its normal unit vector − → u n . At each point of the surface, the 3D electric field is described in the right-handed coordinate system defined by the normal unit vector − → u n , the unit vector in the propagation direction − → u p and in the direction orthogonal to − → u p within the 2D plane − → u ⊥ , see Deno et al. (2017). Since a 2D MEA was used in this study, only the electric field components in the plane defined by − → u p and − → u ⊥ were estimated: traveling wave is assumed to be locally plane and homogeneous within each clique of the MEA. Under that assumption, the relation between the local electric field and the measured b-EGMs is given by the following linear system: where b(t) is a vector containing the available b-EGMs, D is a matrix of interelectrode distances, and ̺ ∈ { , △ } denotes the clique type, square ( ) or triangular ( △ ), respectively. In case of square clique, as the one illustrated in Figure 3A is a 2×6 matrix, reported below. Similarly, in case of triangular configurations, as the clique 1 in Figure 3B matrix, also shown below.
(16) Expressions similar to D △,1 were obtained for the other threeelectrode configurations of the clique shown in Figure 3B. Note that with this formulation, the electric field is taken at the clique center, assuming that there are no delays between b-EGM components.
From Equation (15) we see that the projection of the E-field along the inter-electrode distance results in the corresponding b-EGM. Since measured b-EGMs will be affected by noise and model mismatch errors, we can use the least-squares criterion to estimate the electric fieldÊ ̺ (t) at the center of each clique: The evolution over time of theÊ ̺ (t), estimated at the center c of each ̺-type clique, describes a loop trajectory in the 2D plane.
Under the assumption of locally plane and homogeneous wave, E ⊥ (t) ≈ 0 within the clique, so thatÊ ̺ (t) should have a dominant direction and the loop should lie on a straight line along the propagation direction − → u p (Deno et al., 2017). On the other hand, if the wave is not plane,Ê ̺ (t) will describe a bidimensional loop and propagation cannot be characterized by the projection on a single direction.

Electric Field Estimates at Each Clique
If N×N is the size of the MEA, there will be (N − 1)×(N − 1) = N 2 − 2N + 1 estimates ofÊ(t), in the case of square cliques, and four times this value in the case of triangular configurations. Using Equation (17) and the definitions in Equations (4, 16), the two componentsÊ x (t) andÊ y (t) of the electric field can be expressed as linear combinations of the two b-EGMs along x-and y-directions, respectively: where the approximated terms make use of Equations (7-10).
Equations (18)(19)(20)(21)(22) show that the two electric field components within each clique have a different amplitude, depending on the propagation angle θ . In addition, for triangular cliques, there is also a delay between both components, which is also a function of the propagation angle θ and the conduction velocity magnitude v. Such delay represents the difference in the wavefront arrival time to the middle point of each bipole.

Time Alignment of b-EGMs
In the formulation of Equation (15), no delay is considered among the b-EGMs. However, when estimating the electric fields making use of measured b-EGMs, the different bipolar signals are activated at their corresponding wavefront arrival times, thus having a delay among them, as expressed in Equations (7-12).
For the square cliques, from Equation (18) we found that there are no delays betweenÊ x (t) andÊ y (t), for plane wave propagation and to that level of approximation. Consequently, theÊ (t) loop should lie on a one-dimensional trajectory. However, each componentÊ x (t) andÊ y (t) in Equation (18) is the summation of two terms which have a delay between them: Equations (7-10), respectively. It is well known that the summation of delayed components result in smoothed estimates, (Sörnmo and Laguna, 2005), which can lead to underestimateÊ (t) at the peak of the signal. Therefore, the first-order Taylor series approximation used in this work to arrive to Equation (18) no longer applies. This phenomenon can be observed in the second graph of row in Figure 5D within the Results section, where the E-field estimate is underestimated at the peak of the signal.
For the triangular cliques, Equations (19)(20)(21)(22) show that, in general, there is a delay betweenÊ x (t) andÊ y (t) components. This misalignment will be reflected at theÊ △ (t) loop trajectory, which thus will lie in a bidimensional plane rather than on the unidimensional line in the wave propagation direction. To avoid these phenomena and to have well synchronizedÊ x (t) andÊ y (t), in this study we propose a modified version of the least squares estimator of E ̺ (t), where the b-EMGs in the clique are previously time-aligned. For that purpose, both in square and in triangular configurations, the time delay τ mn between each independent bipolar signal b mn (t), mn ∈ {12, 34, 13, 24} and the b-EGM with the highest amplitude within the clique b max (t) was estimated by maximizing their cross-correlation: Each b-EGM b mn (t) was then advanced by τ mn so that it is aligned with b max (t). Figure 5 shows examples of E-field estimated with and without previous alignment.

Estimation of Propagation Angle
OP-EGM approach estimates propagation direction by finding the angleθ which maximizes the cross-correlation between the first time derivative of the unipolar voltage signal φ ′ (t − t c ), associated to the clique c, and the projection of the electric field in that directionÊ p (t). The unipolar voltage signal φ ′ (t − t c ) at a given clique corresponds to a locally estimated u-EGM, u ′ c (t), eventually delayed by In the latter equation, u ′ c (t) must be estimated for each clique as explained later. In the coordinate system depicted in Figure 1D, the E-field Therefore, the projection of E(t) along the propagation direction is given by: The estimated conduction angleθ and the possible delayτ c will be found by replacing the E-field E(t) with its estimate according to Equations (18-22) and solving: Ifτ c is the estimated delay τ c maximizing the previous crosscorrelation, the propagation angle estimate is: For square cliques, considering the estimates of the E-field components for plane wave propagation in Equation (18), the expressions in Equation (3) and the fact that u ′ Using a similar procedure for triangular cliques, as clique 1 in Figure 3B, considering the fact that u ′ c (t) and u ′′ c (t) are orthogonal and using the approximation in Equation (19), we obtain that the propagation angle estimate when no previous time alignment of b-EGMs is performed can be expressed as: This result points out that the estimate of θ is biased in general for triangular cliques. Nevertheless, if the electric field componentŝ E x (t) andÊ y (t) are estimated from properly aligned bipolar signals, then we obtainθ = θ . For cliques corresponding to triangles 2, 3 and 4, similar expressions are obtained with only changes in the signs of the sin(θ ) and cos(θ ) terms.
In the reference omnipolar method introduced by Deno et al. (2017), it is not specified which of the unipolar signals measured in a clique is used as the u ′ c (t) in Equation (25). In this work, when reproducing the method proposed in literature, we take the signal from lower left corner electrode when using square cliques, so that u c (t) = u 1 (t). For triangular cliques, the signal from the electrode corresponding to the vertex of the right angle is chosen, so that u c (t) = u k (t), k = {1, 2, 3, 4} for each of the four triangles indicated with the same k-th index, respectively, in Figure 3B.
Alternatively, we propose in this work a more robust approach, where the estimation of u ′ c (t) to be included in Equation (25) involves all the u-EGMs within the clique. For the case of square clique: 1. each of the four u-EGMs is aligned with respect to the one with the highest amplitude, denoted as u max (t), by computing the maximum cross-correlation between them. In this way, the four aligned u-EGMs u k (t − τ k ), k ∈ {1, 2, 3, 4}, are obtained; 2. the average signalū(t) is calculated by averaging the four aligned u-EGMs u k (t − τ k ); 3. the time derivative ofū(t) is computed and used as estimate of u ′ c (t) in Equation (25), u ′ c (t) = ∂ū(t)/∂t. For triangular cliques, steps 1. and 2. are performed by just considering the three u-EGMs at the electrodes forming the clique.
Step 3. is the same as for square cliques.

Estimation of Conduction Velocity
According to the omnipolar method, the conduction velocity v in each clique is estimated by comparing the amplitudes of φ ′ (t −t c ) andÊ p (t), as they are expected to be proportional under the assumption of locally plane wave:Ê p (t) = −φ ′ (t − t c )/v, if the b-EGMs within the clique are aligned (see Equations 18-22). The conduction velocity v can be estimated by computing the ratio of the peak-to-peak (pp) amplitudes of φ ′ (t − t c ) andÊ p (t), as proposed by Deno et al. (2017): In order to avoid the large sensitivity to noise that be expected from an estimation based on peak-to-peak amplitudes, we propose here a more robust estimate of v as the ratio of the For square cliques, by using Equation (18) in Equation (24) and if the propagation direction is well estimated (θ = θ ), simple calculations lead to the expected: For triangular cliques (e.g., clique 1), using the approximations with delays in Equation (19) forÊ x (t) andÊ y (t), and assuming that the propagation direction is well estimated (θ = θ ), we have:Ê This equation shows that the misalignment ofÊ y (t) andÊ x (t) components in triangular cliques has an effect in the conduction velocity estimate, even if propagation direction is estimated without error. Nevertheless, if the b-EGMs are time aligned as proposed in this work (i.e., if there are not delays between E y (t) andÊ x (t)), the approximation in Equation (31) reduces to an expression analogous to Equation (30), allowing a proper estimation of the conduction velocity.

Omnipolar Signals for Voltage Maps
In this work, OP-EGM signals ϕ i,j,(q) (t), where q ∈ {1, 2, 3, 4} refers to each of the four possible triangles when using triangular clique configuration, were defined by projectingÊ(t) within the corresponding clique onto the following directions: • The direction − → u me that maximizes the excursion of theÊ(t) loop within the depolarization interval. In that case, the OP-EGM signal ϕ me i,j,(q) (t) along the loop maximal excursion (me) direction was derived as: ϕ me i,j,(q) (t) = d ·Ê(t) · − → u me . It should be noted that the peak-to-peak amplitude of ϕ me i,j,(q) (t) is equivalent to the peak-to-peak amplitude E max of the electric field, proposed by Deno et al. (2017), multiplied by the interelectrode distance within the clique.
• The principal direction − → u ν obtained by the principal component analysis (pca) of the loop described byÊ(t) within the depolarization interval. The resulting OP-EGM signal ϕ pca i,j,(q) (t) was derived as: ϕ pca i,j,(q) (t) = d ·Ê(t) · − → u ν . In this approach, we also obtained the voltage projected in the orthogonal direction − → u ⊥ ν (corresponding to the second principal component): Let us consider ϕ me i,j,(q) (t) under the plane wave assumption. For square cliques, if maximal excursion occurs in the propagation direction and conduction angle is well estimated, we will have: It should be noted that there are no second-order terms in Equation (32), unlike the expressions of b-EGMs (see Equations 7-12). Again, alignment prevents attenuation at theÊ(t) peaks and consequently yields a better estimate. For triangular cliques under the same assumptions (let us consider clique 1, for instance), there are indeed second order terms in ϕ me i,j,1 (t), affecting omnipolar signal amplitude. Nevertheless, expressinĝ E y (t) andÊ x (t) components with the delay approximation in Equation (19) and time aligning them, we obtain the aligned version ϕ me,a i,j,1 (t): whose amplitude is properly estimated. Examples of the omnipolar EGM signals here introduced are observed in Figure 6, for both clique configurations and the three MEA orientations, with and without the previous alignment of b-EGMs.

Voltage and Conduction Velocity Maps
For each OP-EGM signal introduced in this work ϕ s i,j,(q) (t), s ∈ {me, pca, pca⊥}, we computed the peak-to-peak amplitude V o-s i,j,(q) , which was used to build voltage maps. In the pca approach, we also created maps based on the root sum square of the voltage in the first and second principal components, V o-pca−r i,j,(q) : We also built bipolar maps based on the peak-to-peak amplitudes V b-x i,j and V b-y i,j of the b-EGMs in each of the two catheter directions (i.e., b x i,j (t) and b y i,j (t)). The root sum square (r) of both values, as well as their maximum (m) were also considered so as to have a performance reference based on bipolar signals: All the voltage maps were computed for each orientation of the MEA, and by using both clique configurations analyzed in this study. Omnipolar voltage maps were also performed from aligned b-EGMs. We obtained color-coded maps where a pixel corresponds to a clique (or to an electrode pair, in case of bipolar maps). Regarding conduction velocity (v = v − → u v ), we created maps for each catheter orientation, each clique configuration, with and without alignment of b-EGMs. More specifically, each map consists of color-coded v-values and arrows representing angular direction θ of the estimated velocity vectors v in each clique.
Plotted v-values were estimated by using both OP-EGM approach of the Equation (28) proposed by Deno et al. (2017)

and the modified version MOP-EGM introduced in this work in Equation (29).
It should be noted that maps performed with square cliques present the same resolution as bipolar maps, both providing N 2 −2N +1 pixels when processing the whole MEA. On the other hand, since triangular cliques result in a total of 4×(N 2 − 2N + 1) measurement configurations (and consequently, pixels) within the MEA, maps performed with this configuration present higher spatial resolution.

Voltage and Conduction Velocity Maps From Noisy Signals
In order to assess the sensitivity to noise of voltage and conduction velocity maps, simulated b-EGMs were corrupted with noise obtained from real b-EGMs recorded at Hospital Santa Marta, Lisbon. Six hundred different noise segments were extracted from b-EGMs recorded with a PentaRay R catheter (Biosense-Webster, Inc., Diamond Bar, CA, USA) in 27 different mapping points at intervals with no recorded EGMs. All noise segments were normalized in order to guarantee the same power level, coinciding with the observed average power. One hundred different realizations of this recorded noise were randomly added to each of the simulated b-EGMs b x i,j (t) and b y i,j (t) along x and y directions of the MEA. As a result, we generated one hundred noisy bipolar signals along each of the two directions, b x,n i,j (t) and b y,n i,j (t). For each realization n, the different modalities of voltage and conduction velocity maps were computed, obtaining 100 different maps from noisy b-EGMs. These mapping strategies were tested regarding fibrosis detection, as well as voltage map reproducibility. The previous procedure was repeated for six different noise levels, scaling the noise realizations so that they had standard deviations σ n ∈ {3, 6, 14, 28, 42, 55} µV. The standard deviation of the recorded noise ranged from 1.2 µV to 11 µV, with an average value of 2.4 µV, lower than those reported in Unger et al. (2019). We simulated noise levels ranging up to 55 µV which approximately corresponds to the 95th percentile of the bipolar noise levels reported in Unger et al. (2019), in order to test the performance of the mapping strategies at challenging situations. An example of noisy b-EGM b y 14,2 (t) with noise level σ n = 14 µV is shown in Supplementary Figure 2C.

Maps Performance Evaluation for Fibrosis Detection
In order to quantitatively evaluate the ability of the different mapping strategies to identify the fibrotic area (i.e., to discriminate pixels within the fibrotic patch from those in the normal tissue), receiver operating characteristic (ROC) curves were used. Two ground-truth masks of fibrosis were created for bipolar and omnipolar maps, having the resolution of square and triangular cliques and indicating whether a clique lay within the fibrotic or the normal tissue. A 14 × 14 ground-truth mask was used to evaluate bipolar and omnipolar maps with square cliques. A value of 1 was assigned if the four electrodes defining the square clique lay in the fibrotic area, and a value of 0 was assigned if all of them lay in the non-fibrotic area. Cliques with some electrodes inside and some outside the fibrotic patch were not considered in the evaluation. Similarly, a 28 × 28 ground-truth mask was used to assess omnipolar maps performed with triangular cliques, considering if the three electrodes defining the triangular clique lay in fibrotic or non-fibrotic tissue.
Both ground-truth masks used in this study are illustrated in Figures 4A,B. For each voltage and velocity map, ROC curves relating sensitivity and specificity were obtained by varying the threshold for fibrosis identification (Metz, 1978). In this analysis, 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 accuracy (ACC) was defined as the number of cliques correctly identified as fibrotic or non-fibrotic divided by the total number of cliques. The maximum ACC was used as a measure of the overall fibrosis detection ability of a given map. In addition, the sensitivity and specificity with the threshold achieving the maximum ACC were also computed.

Maps Performance Evaluation for Reproducibility of the Reference Voltage
The ability of the voltage mapping strategies in reproducing the spatial distribution of reference maps was also evaluated. We created peak-to-peak amplitude maps of u-EGMs, and resampled them to match the spatial resolutions of bipolar and omnipolar pixel maps (with both square and triangular cliques). To quantify the agreement of each voltage map with the substrate characterization given by the reference map, we used both Pearson's ρ p and Spearman's rank ρ s correlation coefficients. A p ≤ 0.05 was required to have statistical significance. Reference voltage maps were obtained by following two sequential steps: 1. A bicubic interpolation of the N×N peak-to-peak amplitude map of u-EGMs u i,j (t), (i, j) ∈ {1, ..., N}, was performed FIGURE 4 | (A) 14×14 and (B) 28×28 ground-truth masks for fibrosis detection, for square and triangular cliques, respectively. (C) 15×15 peak-to-peak amplitudes map from u-EGMs and (D) 57×57 map resulting from bi-cubic interpolation of (C). 14×14 (E) and 57 × 57 (F) amplitude maps for evaluating voltage reproducibility of bipolar and omnipolar maps.
placing three interpolated values between any two adjacent electrodes, along both MEA directions, thus increasing the total number of pixels to (N + 3 · (N − 1)) × (N + 3 · (N − 1)). The original and interpolated maps for N = 15 are shown in Figures 4C,D, respectively, when MEA is oriented with = 0 • . 2. The interpolated pixel map was later resampled so as to obtain (N − 1) × (N − 1) and 2 · (N − 1) × 2 · (N − 1) reference maps, matching the dimensions of bipolar and omnipolar maps with square cliques, and omnipolar maps with triangular cliques, respectively. Resulting resampled pixel maps from N = 15, for the case when = 0 • , are illustrated in Figures 4E,F.

Error Evaluation for Noisy Maps
Additionally, in order to evaluate the effect of noise on voltage and velocity maps, we computed the root mean square error (RMSE) between each noisy map and the respective noise-free map.
In case of voltage maps, we calculated: where s ∈ {b-m; b-r; o-me; o-me, a; o-pca-r; o-pca-r, a}. Analogously, for velocity maps, we computed: with s ∈ {o; o-m; o-m, a}.
We also studied the effect of noise on the estimation of the propagation direction with both the original θ o and the modified θ o-m omnipolar method. Angle differences, ǫ θ , between propagation direction maps obtained from noisy realizations and those obtained from noise-free signals (considered as the reference) were computed for each omnipolar approach and noise level. The mean and standard deviation of the angle differences in each map were obtained and averaged for all noise realizations.

RESULTS
The effect of the proposed alignment is evident in both upper and lower panels. In square cliques, the E-field peak value for the aligned versions is clearly higher than for the non-aligned ones. On the other hand, in triangular cliques the alignment of b-EGMs proves to improve the estimation of the propagation direction. These is clearly seen especially in those areas where the local propagation direction differs from the dominant directions of the MEA [see panels (B), (C), (D) and (E)]. Figure 6 shows the OP-EGM signals estimated with the different proposed approaches at the same electrode positions(i, j) as in Figure 5: ϕ s i,j,(q) (t), where s ∈ {me, pca, pca⊥}, in square and triangular cliques, with and without alignment of b-EGMs. These results show that the OP-EGM signals produced by the two modalities are quite similar in most of the cases. Square cliques reveal that ϕ me i,j (t) is equivalent to ϕ pca i,j (t), both when the b-EGMs are aligned and when they are not. In triangular cliques, this is also true in areas where the propagation direction is parallel to one of the main directions of the MEA.
The different bipolar and omnipolar voltage mapping strategies are illustrated in Figure 7 for = 0 • , without noise (left columns) and with a noise level of σ n = 28µV (right columns). In the latter case, only one of the 100 noisy realizations was shown. The omnipolar maps are presented without and with the alignment of b-EGMs. These results reveal that combined bipolar maps obtained as the root sum square or the maximum voltage of both MEA directions, V b-r and V b-m , shown at row A, present better fibrosis detection performance than omnipolar maps without alignment, illustrated at upper and lower rows in Figure 7B for the two approaches me and pca, respectively. However, omnipolar maps with previous time alignment of b-EGMs (upper and lower rows in Figure 7C for the same approaches as in Figure 7B) show results comparable to combined bipolar maps. Omnipolar mapping strategies also present greater correlation with the reference map when time alignment is applied. When b-EGMs are affected by noise, omnipolar voltage maps are more robust than bipolar maps, improving their performance in fibrosis discrimination and correlation with their reference when previous alignment is performed.
Velocity maps obtained with the reference OP-EGM method (row A) and with the proposed approach (without (row B) and with (row C) alignment of b-EGMs) are presented in Figure 8, at the same MEA orientation and with the same noise level as in Figure 7. These maps reveal better performance of the MOP-EGM method (CV o-m and CV o-m,a ) than the standard approach (CV o ). Moreover, the standard approach tends to overestimate conduction velocity values, when compared with the MOP-EGM method.
AUC, ACC, Pearson ρ p and Spearman ρ s correlation values of the proposed mapping strategies for noise-free signals are shown in Table 1, together with the threshold value with maximum detection accuracy for each map modality. Note that performance was computed by considering jointly the three MEA orientations with respect to the propagation direction. Bipolar voltage maps V b-m and V b-r identify the fibrotic area with an ACC of 96 and 95%, respectively, whereas the OP-EGM based voltage maps reach a maximum ACC value of 93%, being consistent with results in Figure 7. Regarding correlation of the voltage maps with the reference maps, both omnipolar strategies (me and pca) achieve values ρ p = 0.87 and ρ s = 0.87 when square cliques are used and b-EGMs alignment is performed. As a reference, V b-r and V b-m present correlations of ρ p = 0.86 and ρ p = 0.83, respectively. Regarding CV maps, the MOP-EGM version proposed in this work (with b-EGMs alignment) presents indeed comparable performance in fibrosis detection to combined bipolar voltage maps (e.g., ACC = 96% for CV o-m,a with square cliques), and better performance than both omnipolar voltage maps and than the original OP-EGM approach to estimate conduction velocity (CV o ).
Similar results are observed in Figure 9, which shows how fibrosis detection (in terms of the ACC value, Figure 9A) and voltage map fidelity (in terms of Pearson's correlation coefficient ρ p , Figure 9B) are affected when b-EGMs are corrupted with different noise levels, jointly considering the three MEA orientations. As observed in panel A, omnipolar voltage and CV maps are less affected by noise than bipolar maps for noise levels greater than or equal to σ n = 28µV. Note that the proposed time alignment of b-EGMs improves fibrosis detection  , are the most robust in reproducing the reference unipolar voltage maps in the presence of noise. They achieve Pearson correlation coefficients >0.80 for all noise levels and present smaller interquartile ranges than bipolar maps.
The RMSE between maps computed from noisy b-EGMS and those computed with noise-free b-EGMs are depicted in Figures 10A,B, respectively, for each noise level. As in previous results, all three MEA orientations are considered jointly. The RMSE values for omnipolar voltage maps are all lower than 0.2 mV, whereas they achieve 0.3 mV in bipolar maps for the highest noise level. As for the CV maps, the MOP-EGM method shows lower RMSE values than the standard OP-EGM approach, in agreement with what was found in Figure 8. Table 2 reports the mean and standard deviation of the error caused by each noise level in the propagation angles estimated with both omnipolar approaches θ o and θ o-m for square and triangular cliques, considering jointly all three MEA orientations. These results show that propagation direction maps performed with square cliques are less affected by noise than those performed with triangular cliques, regardless of the omnipolar approach. However, for noise levels lower than σ n = 28µV, direction maps with the MOP-EGM approach θ o-m exhibit a smaller mean bias than those using the original method θ o .
Finally, Figure 11 shows the effect of the MEA orientation on omnipolar voltage and velocity maps (V o-me , V o-pca-r , V o-me,a , V o-pca-r,a , CV o-m and CV o-m,a ) performed with square cliques.

DISCUSSION
In this work, we investigated the performance of different substrate and velocity mapping modalities based on the OP-EGM method, and on the here proposed MOP-EGM approach, in characterizing and detecting fibrotic areas. For that purpose, we used a 2D multi-electrode array over an atrial tissue simulated with the Courtemanche cellular model.
The OP-EGM methodology was presented in Deno et al. (2017) as a novel approach allowing characterization of myocardial substrate and propagation pattern, making possible real-time high-density maps less sensitive to catheter orientation than current bipolar strategies. It showed to solve complex collision and fractionation wavefront patterns in animal subjects, thus providing coherent voltage maps even during atrial fibrillation (Haldar et al., 2017). The OP-EGM approach was proposed to overcome the well-known weaknesses affecting   bipolar voltage mapping, which is the cornerstone for the identification of fibrotic atrial substrate, and consequently for low-voltage-guided catheter ablation strategies in AF (Haldar et al., 2017). The efficacy of omnipolar EGMs was also proved in ex vivo  and in vivo (Porta-Sánchez et al., 2019) ventricular mapping, for delineation of both healthy and infarcted areas. OP-EGM approach-based mapping strategies represent an alternative to the more timeconsuming local activation times mapping, allowing beat-to-beat indication of the wavefront propagation direction (Deno et al., 2020). However, we observed in our work that OP-EGM still presents some orientation-dependent behavior affecting fibrosis detection, voltage and CV estimations, and we propose an alternative, referred as MOP-EGM method, partially overcoming these problems. We started by analyzing maps computed from noise-free b-EGMs. Our results show that the previous alignment of the b-EGMs in each clique (as proposed in the MOP-EGM method) always improves the fibrosis detection ability of the omnipolar voltage and conduction velocity maps. The latter provide a characterization of the atrial substrate comparable to bipolar voltage maps combining b-EGMs along x and y directions. On the other hand, bipolar maps V b-x and V b-y in each of the MEA main directions revealed ACC values of 69 and 91%, respectively. Moreover, CV maps computed with the MOP-EGM approach present better performance than both omnipolar voltage maps Frontiers in Physiology | www.frontiersin.org and CV maps computed with the original OP-EGM method. It must also be pointed out that, using the omnipolar approaches, conduction velocities are estimated without the need of detecting local activations. Therefore, the sensitivity to inaccuracies in the determination of local activation times is avoided. Most studies to detect fibrosis rely on bipolar voltage mapping (Rodríguez-Mañero et al., 2018;Yamaguchi et al., 2019). Other works demonstrate a high correlation with respect to the use of u-EGMs based maps, suggesting no relevant clinical impact of using b-EGMs over u-EMGs (Nairn et al., 2020). Future comparison of unipolar maps with the here proposed MOP-EGM approach will elucidate to what extend that statement remains valid.
The present work also aims at studying the sensitivity to noise of the different mapping approaches. Our observations suggest that the voltage and conduction velocity maps performed from aligned b-EGMs are more robust against noise than both bipolar maps and their unaligned versions, also showing better performance in discriminating fibrosis and in reproducing voltage maps for the same noise level.
In the present study we also addressed the question of the sensitivity of the OP-EGM method to the MEA orientation. Although the fibrotic patch was well discriminated from the healthy tissue by all omnipolar mapping strategies considered in this study, both voltage and velocity OP-EGM based estimates were affected by the relative orientation between MEA and tissue fibers. Regarding voltage maps, when the MEA is not oriented in the same direction as tissue fibers (as when = 0 • ), we observe a reduction of the omnipolar voltage values. In addition, in tissue areas where the wavefront propagation is curved, the electric field is overestimated when the local propagation direction becomes parallel to one of the two main directions of the catheter (as when = 45 • ). As this happens in an area of healthy tissue, it causes a higher overall accuracy in fibrotic patch detection. The reverse occurs when these curved areas are oblique to the MEA (as in = 0 • ), providing underestimated voltages and decreased fibrotic detection accuracy.
As for the orientation dependence of velocity maps, the opposite behavior is observed, with velocity estimation increasing for = 0 • . This is explained by the fact that v, in Equation (29), has an inverse relationship with the magnitude of the estimated E-field, while in voltage estimation (equations in subsection 2.9), the dependence is direct. In addition, in those areas where the wave propagation is curved, velocity is overestimated when the MEA is in the same direction as the tissue ( = 0 • ). Therefore, the fibrotic patch is globally better detected for this configuration. On the other hand, if catheter becomes oblique to the tissue ( = 0 • ), areas of curved propagation appear underestimated, thus reducing the overall fibrosis detection performance for those configurations. When b-EGMs are time-aligned, both MOP-EGM based voltage and velocity maps continue to be sensitive to the orientation between MEA and tissue, but to a lesser extent if compared to their unaligned versions. As a result, improved fibrosis detection, better voltage estimations and less degradation by noise are observed in average when using the MOP-EGM approach. The remaining dependence with is not predicted by the mathematical model based on the plane-wave assumption, but the sampling frequency of the b-EGMs (1 kHz) could be behind this observation.

LIMITATIONS
The basic assumption of plane and homogeneous propagation within each clique represents an intrinsic limitation of the omnipolar methodology, requiring the use of high spatial resolution multi-electrode catheters (like the one considered in this study). Moreover, the hypothesis of a plane myocardial surface in the area defined by the clique is questionable in non-compacted myocardial regions or in the atrial appendage. Nevertheless, those areas will be identified with 3D electroanatomic mapping systems (Deno et al., 2017).
In this work, we simulated a very simple propagation pattern in a 2D atrial model, which does not perfectly reproduce the real 3D situation. Therefore, the results here obtained need to be confronted in future studies with those from other simulation configurations (such as non-plane or multiple wavefronts in atrial fibrillation, as well as patchy fibrosis models). Moreover, the proposed MOP-EGM maps need to be evaluated over real data, where the effect of electrode size, inter-electrode distance and variable tissue-electrode contact will be present. The validation of these methods in-vivo is, though, challenging, since late gadolinium-enhanced magnetic resonance imaging represents the only non invasive available reference for atrial fibrosis today, and even their utility is under debate (Caixal et al., 2020). In-vivo/Ex-vivo animal experiments represent an alternative for future works aiming to establish if the simulation-proved superior performance of MOP-EGM translates to clinical counterparts.

CONCLUSION
In this work we showed that voltage mapping strategies based on the MOP-EGM method are able to discriminate fibrotic from healthy tissue. For low noise levels, they attain comparable performance to maps obtained by combining the b-EGMs amplitudes along the two directions of the catheter, but clearly outperform bipolar maps when noise is added. On the other hand, omnipolar CV maps reveal higher fibrosis detection accuracy than voltage maps, specially when the modifications proposed in this work were applied. The fibrosis detection performance of voltage and velocity maps benefit from prior intra-clique time alignment of b-EGMs as proposed in the MOP-EGM approach. In addition, the use of square cliques outperforms the triangular ones.
As for the different approaches to estimate the omnipolar voltage from the E-fiel loop, both the one based on maximal excursion (me) and the one based on principal component analysis (pca) yield similar results, in terms of detection accuracy and correlation with unipolar voltage maps. Similar conclusions are drawn when b-EGMs are affected by noise extracted from real signals.
Both OP-EGM and MOP-EGM strategies considered to estimate the propagation angle revealed smaller sensitivity to noise than methods based on bipolar voltage maps. The use of square cliques results in lower variance than the use of triangular ones when estimating the propagation angles in noisy situations.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Central Lisbon University Hospital Centre (CHULC) Ethics Committee. The patients/participants provided their written informed consent to participate in the study.