Improved local activation time annotation of fractionated atrial electrograms for atrial mapping

Background: Local activation time (LAT) annotation in unipolar electrograms is complicated by interference from nonlocal atrial activities of neighboring tissue. This happens due to the spatial blurring that is inherent to electrogram recordings. In this study, we aim to exploit multi-electrode electrogram recordings to amplify the local activity in each electrogram and subsequently improve the annotation of LATs. Methods: An electrogram array can be modeled as a spatial convolution of per cell transmembrane currents with an appropriate distance kernel, which depends on the cells’ distances to the electrodes. By deconvolving the effect of the distance kernel from the electrogram array, we undo the blurring and estimate the underlying transmembrane currents as our desired local activities. However, deconvolution problems are typically highly ill-posed and result in unstable solutions. To overcome this issue, we propose to use a regularization term that exploits the sparsity of the first-order time derivative of the transmembrane currents. Results: We perform experiments on simulated two-dimensional tissues, as well as clinically recorded electro- grams during paroxysmal atrial fibrillation. The results show that the proposed approach for deconvolution can improve the annotation of the true LAT in the electrograms. We also discuss, in summary, the required electrode array specifications for an appropriate recording and subsequent deconvolution. Conclusion: By ignoring small but local deflections, algorithms based on steepest descent are prone to generate smoother activation maps. However, by exploiting multi-electrode recordings, we can efficiently amplify small but local deflections and reveal new details in the activation maps that were previously missed.


Introduction
Atrial electrograms (EGMs) recorded during high resolution atrial mapping as well as their corresponding activation maps, facilitate the identification and localization of potential triggers and substrates perpetuating atrial fibrillation (AF). They can also be used to guide cardiologists through treatment approaches such as EGM based ablation therapies [1].
An EGM recording is an aggregation of the per cell activities that are in the neighborhood of an electrode. How strong each individual cell contributes to the total aggregated potential depends on the distance of the cell to the electrode, as well as on the electrode's dimension [2]. The depolarization wavefront propagation in a rather homogeneous tissue with a uniform excitation wavefront results in a smooth stereotype atrial activity: a positive spike followed by a negative deflection. In these cases, the local activation time (LAT) of the cells that are right under the electrode coincides with the maximum negative deflection or the steepest descent (SD) of the electrogram. * Corresponding author. However, if the tissue is inhomogeneous or multiple excitation wavefronts are propagating in the electrode's neighborhood, the recorded EGM may contain nonlocal (far-field) deflections from the so called ''distant active membrane'' [3]. In such cases, the steepest descent annotated as the LAT might belong to a far-field excitation and not to the local activity. This makes the estimation of LATs prone to errors and will negatively affect the estimation of other parameters that depend on LATs, such as conduction velocity or tissue conductivity [4]. This, in turn, complicates the analysis of EGM recordings and the underlying atrial substrate for AF by adding confusion about the origin of observed components [5]. As shown in a study on bipolar electrograms in [6], 67% of complex fractionated electrograms represent nonlocal activities and are thus not consistent with spatial disorganization in the tissue.
To attenuate the effect of nonlocal activities, bipolar electrograms are alternatively used in many studies. However, these recordings suffer from serious drawbacks including sensitivity to propagation direction and providing no useful information about the potential distribution in the electrogram. Many other studies directly focus on improving LAT estimation in fractionated atrial electrograms. A review of such methods can be found in [7]. Some of these approaches try to make use of spatial information offered by multi-electrode recordings. These include application of cross-correlation, spatial filtering, spatial gradient or spatial deconvolution, among which only the spatial deconvolution is formulated according to an electrophysiological model of the tissue. It benefits from the fact that a multi-electrode electrogram recording can be modeled as a convolution of the per cell transmembrane current with a spatial filter that depends on the electrode's diameter and its distance to the cells [8].
In this paper, we also exploit the convolutive electrogram model for a better annotation of the local deflection (and respectively the local activation time) in fractionated electrograms. We first present a model for the EGM array that reflects the spatial convolution. By employing this model, we formulate a deconvolution problem to estimate the local activities. However, the deconvolution problem is under-determined and results in unstable solutions. To overcome this issue, some studies employ prior knowledge such as planar wavefront propagation [8].
Other studies use spatial interpolation and filtering to artificially introduce more data points [2]. None of the assumptions in these approaches are satisfactory from a physiological point of view. In this paper, we propose to use a regularization based on sparsity of the first-order time derivative of the electrograms to solve the inverse problem of local activity estimation. We solve the proposed deconvolution problem using the Split Bregman method [9]. Overall, the approach is computationally efficient, does not suffer from boundary artifacts, and can also perform with incomplete observations. Our results on simulated and clinically recorded data show that the proposed approach can efficiently amplify the local activities in an EGM and outperforms the reference approaches in estimation of the LAT.

Action potential propagation
An electrogram is the record of changes in the potentials of the many cells surrounding an electrode. These changes in per cell potentials are the result of action potential propagation in the tissue. This propagation in a two-dimensional (2D) tissue can be modeled using the following reaction-diffusion equation [10] ( , , ) where ( , , ) is the per cell potential at location ( , ) and time , = 1 μF cm −2 is the total membrane capacitance, st is the stimulus current, and ion is the total ionic current computed according to the Courtemanche model in [11]. ( , , ) = −1 ∇⋅ ( , )∇ ( , , ) is the transmembrane current, with = 0.24 μm −1 the cellular surfaceto-volume ratio, and ( , ) the intracellular conductivity tensor.

Electrogram model
An electrogram can modeled as the weighted summation of per cell transmembrane currents. The weights are equal to the inverse of the cell-to-electrode distance. The distance between a cell at position ( , ) and an electrode at position ( , ) and a (constant) height 0 above the 2D tissue equals √ ( − ) 2 + ( − ) 2 + 2 0 . The electrogram can thus be modeled as [10] ( , , ) = for = 1, 2, … , , where is the total number of electrodes,  is the area in which the modeled cells are located, ( , ) is the area variable, and is the constant extra-cellular conductivity. Eq. (2) represents a 2D spatial convolution of ( , , ) with a distance kernel 0 given by A plot of 0 is shown in Fig. 9(a). The spatial sampling implemented by the electrode array is modeled by a sampling operator 0 , where ( ) is a Dirac delta impulse. The measurement model in Eq. (2) can thus be written as where * * denotes the 2D spatial convolution. Note that we have assumed the electrode diameter can be neglected in the model for 0 . We will motivate this assumption in Section 3.6.
To be able to invert the model, we discretized Eq. (4) in space as well as in time. In space, we use source clamping and replace each block of cells in the three dimensional tissue with a modeled ''cell'' on a uniform 2D grid of cells with cell-to-cell distance and = × the total number of cells. and are the number of rows and columns of the grid, respectively. We also sample in time with sample period and total number of time-domain samples . The discretized convolutive model of the electrogram then becomes where , , are integers that index the sample grid, and = 2 ∕4 contains all constants and will be omitted for simplification. The sampled distance kernel 0 [ , ] is represented by a limited support of size (2 + 1) × (2 + 1) elements. To implement the convolution, the modeled 2D grid of cells need to be extended by cells in each direction. The sampling operator 0 [ , ] selects only the spatial locations on the grid on which we have measurements, and replaces the other locations with zero. This discretized model is easily translated into a matrix model [4].

Transmembrane current estimation
Our aim, in this section, is to estimate the transmembrane currents [ , , ] in Eq. (5), which is both a deconvolution and an interpolation problem. Basically, the deconvolution can be performed by using a loss function that minimizes the least square error between the target value and the estimated value 0 ( 0 * * ) (first term in Eq. (6), below). However, since the number of available electrograms is less than the number of modeled cells and the distance kernel has a lowpass filtering effect, the inverse problem is highly ill-posed and results in unstable solutions. To stabilize the solution, we need to introduce regularization constraints that rule out physically unlikely solutions by employing prior knowledge of the expected solution.
A classical regularization technique in this context is Tikhonov regularization [12], however, it provides only a spatial constraint, and the underlying assumptions on the prior are not strong enough to provide for much interpolation. A priori space-time information such as a specific wave pattern propagation is too specific and does not hold if the EGMs are fractionated. On the other hand, a priori information should be simple enough to allow an efficient problem solving. As shown in many studies, the sharp deflections in an electrogram are of high importance in the analysis of wave propagation and AF.  are its fast temporal deflections. Note that the assumption may not hold for electrograms recorded at areas with continuous electrical activity and no distinct deflection. However, in most types of fractionated electrograms, specially those due to far-field atrial activities, lines of blocks and wave collisions, the number of deflections is still small and the firstorder time derivative of the transmembrane current can be considered sparse. In an optimization framework, this is implemented by imposing an 1 -norm constraint as the regularization function (the second term in Eq. (6), below). The 1 -norm which is the sum of absolute values, is known to induce sparsity in the solution. The resulting regularized optimization problem that we propose is thus given by and is a regularization parameter that is a weight on the importance of the regularization, which will be discussed further in Section 3.7. Due to the coupling between the 1 -norm and the 2 -norm terms, Eq. (6) is a hard problem to solve. An efficient numerical approach to solve 1 -regularized problems is the Split Bregman algorithm [9] which splits the 1 -norm and 2 -norm components by introducing new variables. To efficiently employ the Split Bregman algorithm on this problem, we propose to use two new splitting variables, 1 = 0 * * and 2 = ′ . The new optimization problem will then be where 1 and 2 are the Bregman iterative parameters and 1 and 2 are the penalty parameters. Due to the decoupling, we can now break the problem in Eq. (7) into five easy steps. Each step updates the value of one of the unknown parameters , 1 , 2 , 1 , and 2 . A proper representation of Eq. (7) for structured matrix computations and a detailed description of the five steps taken to solve it, can be found in our previous work in [13], where in a preliminary study, we focused on the proper formulation of the present problem and the derivation of efficient solutions for its steps. In total, the algorithm has a fast convergence rate to a reasonable precision in practice by avoiding costly matrix inversion and performing the computations in the Fourier domain.
Note that the modeled cell size that we use in the inverse problem to estimate the transmembrane currents is much larger than a real cell size and, as mentioned before, each modeled cell in our simulation represents a group of cells in the 3D real tissue. Therefore, the estimated transmembrane currents would actually be more localized electrograms than the exact per cell transmembrane currents. However, in this paper, we will keep referring to them as transmembrane currents.

Strategies for generation of fractionated electrograms
To perform an instrumental evaluation of our proposed approach, we need to generate simulated fractionated electrograms. We use various patterns of heterogeneity for the tissue's conductivity map, based on the literature, to produce fractionation. This assures that the patterns are sufficient to obtain fractionation and the resulting electrograms are representative of real clinical data. The size of each map is 101 × 83 cells, with = 0.6 mm. The simulated patterns, in order of increasing level of heterogeneity, are: 1. A small conduction block: A homogeneous tissue except for two small areas with smaller conduction in the center (see T1 in Fig. 1). Notice that T1 is not representative of real tissue, but serves as a simple example for visualization purposes. 2. Zones of no conduction: In this model the heterogeneities in conductivity are incorporated as a set of randomly positioned lines of blocks, disconnecting the coupling between the cells on the grid [14] (see T2 in Fig. 1). 3. Percolation: the conduction disturbance in this approach is modeled by randomly disconnecting the coupling between some modeled cells and their neighbors [15] (see T3 in Fig. 1). 4. Zones of no conduction and percolation: both approaches are used simultaneously to generate new pattern of heterogeneity (see T4 in Fig. 1).
To model action potential propagation in the simulated tissues, Eq. (1) is discretized and solved using a finite difference method with no flux boundary condition. A more detailed description of the simulation steps and parameters can be found in [4]. The resulting activation maps are shown in the second row of Fig. 1. Each pixel in the activation map represents the true activation time of its corresponding cell which is annotated as the time when the cell's potential reaches a threshold value of −40 mV in the depolarization phase of its action potential. The white pixels belong to the cells that were positioned on a block and did not get activated. Finally, Eq. (2) is used to compute the electrograms recorded by an assumed 15 × 9 electrode array (M=135), with an interelectrode distance of = 3 , covering the 43 × 25 central cells and a constant value of 0 = 1 mm (which will be discussed in Section 3.4). We used the largest kernel size possible, (2 + 1) × (2 + 1) = 83 × 83 to compute the simulated electrograms. The last row of Fig. 1 shows the simulated electrograms recorded at three locations on the tissue denoted by the * . The true LAT of each electrogram which equals the activation time of the cell that is exactly under the electrode, is also denoted by red vertical dashed lines. Due to the heterogeneity in the tissues, the recorded electrograms have multiple deflections. Moreover, as we move from T1 to T4, the level of fractionation increases and the deflections get less distinct.

Transmembrane currents
To visualize the performance of our proposed approach, we start with examples of single electrograms. Fig. 2 shows four examples of fractionated electrograms selected from T1 to T4, respectively. As can be seen in the zoomed activation maps of the 21 × 21 cells surrounding the fractionated electrogram, the slow conductions and blocks in the wave propagation induce fractionation in the recorded electrogram. The first-order time derivatives of the electrograms are plotted in the second row. As can be seen, the local deflection that coincides with the true LAT, is not the steepest descent. These examples show that, although in most studies the steepest descent is used as the LAT, it may not provide reliable results in fractionated electrograms.
We used the 15 × 9 electrogram array and Eq. (7) to estimate the transmembrane current [ , , ]. For simplification the limited kernel size of (2 + 1) × (2 + 1) = 13 × 13 was used in the inverse problem. Since the focus of this study is on introducing the approach itself and not on optimal tuning of its parameters, we used values that yielded good results in LAT estimation in all simulations, which are, = 7×10 −4 , 1 = 1, 2 = 1, and the number of iteration N itr = 100. The custom written MATLAB codes for simulation of atrial electrograms and estimation of LATs using the proposed approach are available in the linked Research Data.
B. Abdi et al. Fig. 3. Two example electrograms in which both SD nd SD fail to estimate the true LAT. The first column shows the zoomed activation maps around the selected electrograms from T2 and T4. The second column shows the electrograms with the true and estimated LAT using both approaches.
The first-order time derivative of the estimated transmembrane currents are plotted in the third row of Fig. 2. As can be seen, the local deflections are amplified and the far-field activities are attenuated. Now as expected, the steepest descent of the transmembrane currents coincides with the true LAT. Fig. 3 also shows two examples where both methods fail to estimate the correct LAT. This highly depends on the wave propagation pattern but mostly happens around areas with a complex propagation pattern, with a strong far-field activity that is very close to the electrode location.
As can be seen, the estimated transmembrane currents do not look as sharp as expected and still some level of fractionation is observed. This is because the estimation of the transmembrane current in the inverse problem is highly ill-posed and there is no guarantee that the resulting electrograms should look like as expected specially in case of fractionation. This can be imposed on the problem by using different regularizations such as imposing sparsity constraint on the transmembrane current itself to provide sharper results or minimizing the error between the estimated currents and a stereotype transmembrane current. However, our goal in this paper is to improve LAT estimation, and our experiments with different regularization terms shows that the current problem formulation in Eq. (6) provides the best results. That is because imposing the sparsity constraint on the first time-derivative of the transmembrane current, helps to preserve the sharp deflections of the transmembrane current (that are later annotated as LATs) while performing deconvolution.

LAT Estimation
Since the transmembrane currents are less affected by the far-field activities, they should provide a better estimation of LATs than the electrograms. To evaluate the performance of our proposed approach (SD ) in LAT estimation, we have compared its results with two reference approaches: (i) steepest descent of the electrograms, and (SD ) (ii) maximum of the spatial gradient of the electrogram array (SG ) [16]. From each of the three conductivity patterns demonstrated as T2, T3 and T4 in Fig. 1, 5 randomly generated tissues were simulated. This Table 1 RMSE (ms) in LATs estimation of the simulated electrograms using the proposed approach SD , the steepest descent SD , and the maximum spatial gradient SG . In the estimation of the error, electrograms that were positioned on the cells that did not get activated were excluded. This resulted in 658, 650, and 649 simulated electrograms in tissue type T2, T3 and T4, respectively. The resulting root-mean-square errors (RMSE) between the true LAT and the LAT estimated by each approach per tissue type are shown in Table 1. On average, in all tissue types, the proposed approach outperforms the reference approaches in LAT estimation.
For a better comparison, Fig. 4 shows the histogram of the absolute errors (larger than 2 ms) in LAT estimation by SD and SD . As can be seen, the number of errors with large values are much larger with SD than with SD . Since, these large errors belong to areas with blocks in the tissue, this shows that a large number of the blocks and heterogeneities are ignored by SD .
As an example, Fig. 5 demonstrates the true activation map of the 15 × 9 electrogram array recorded on T1, as well as the three estimated activation maps using SD , SD , and SG . As can be seen, there is an area with slower conduction in the center of the tissue which causes fractionation in the resulting electrograms. Since this area is small and has a low conductivity, its local deflections are smaller than the nonlocal deflections caused by far-field activities of the neighboring cell. As a result, SD and SG both miss the local deflections, and annotate the far-field activities as the LAT. On the other hand, the proposed approach first amplifies the local activities and then annotates them as the LAT, which helps to preserve the block in the center.   6 also provides another example of a simulated tissue with two extensive lines of block that can potentially lead to a reentry circuit (denoted by the white arrows) in larger scales. The white rectangle in the middle of the tissue denotes the location of the 15 × 9 electrode array. The activation maps estimated using the steepest descent SD and our propose approach SD of the simulated electrograms are also shown. As can be seen, no evidence of the block lines is visible in the SD . This is because the local deflection cased by the wave propagating in between the two lines is very small compared to the strong activity of the neighboring cells. On the other hand, our propose approach SD accurately annotates the small local deflections except for the two electrodes at the boundaries due to the incomplete observations in their neighborhoods.

Electrode height z 0
Throughout this paper, we used a constant value 0 =1 mm for the electrode height, as suggested in [17,18]. In this section, we aim to examine if this value is applicable in our simulation setup and how it affects the simulated electrograms. To do so, we have designed an experiment as follows: (i) first, an 8 × 8 subsection of the clinical electrogram array with equal number of rows and columns is selected. To avoid further complications, a subsection of the array with non-fractionated electrograms is selected. (ii) The activation map of the selected subsection is estimated and interpolated (using cubic interpolation). The resolution is increased by 3 times in both the and the direction. The higher resolution activation map is used as the activation map of the modeled cells whose sizes are one-third of the inter-electrode distance. (iii) To regenerate electrograms based on the higher resolution activation map, the simplified electrogram model in [4] is used. It employs the simplifying assumption that once activated, all cells produce the same stereotype action potential 0 ( ). The stereotype 0 ( ) is estimated from the Courtemanche model for a single cell. The action potential produced by each modeled cell will therefore be ( ) = 0 ( − ), where is the activation time of the cell. (iv) The transmembrane currents are then calculated using ( , , ) = −1 ∇ ⋅ ∇ ( , , ). We further normalize the amplitudes of regenerated and clinical electrograms and ignore all the constants. (v) Finally, Eq. (2) is used to calculate the regenerated electrograms from the transmembrane currents. Since the LATs and other parameters are estimated from the clinical electrograms, the only left parameter that can affect the morphology of the regenerated electrograms is 0 .
As an example, Fig. 7(a) shows the high resolution activation map estimated from the selected 8 × 8 subsection of a clinical electrogram array. Fig. 7(b, c and d) show a clinical electrogram as well as its corresponding regenerated electrograms (using the above approach)   6. Simulated tissue with two parallel vertical lines of block close to the center of the tissue that can potentially lead to a reentry circuit (denoted by white arrows). The white rectangle in the middle of the tissue denotes the location of the 15 × 9 electrode array. The activation maps estimated using the steepest descent SD and our propose approach SD are also shown.
with three different values of electrode height: 0 = 1 mm, 0 = 0.2 mm, and 0 = 5 mm, respectively. The electrogram location is denoted by the red * on the activation map. As can be seen, the effect of 0 is more evident on the steepness of the main deflection and the regenerated electrogram with 0 = 1 mm has the best matching with the clinical electrogram. Although an extension of the above approach can be used for the systematic estimation of 0 , the accuracy of the result will largely depend on the accurate estimation of other underlying parameters including activation map, conductivity map and the stereotype 0 . This makes the correct estimation of 0 in a clinical setup prone to error and almost impractical. Fig. 8 shows two examples of using wrong values for 0 in our proposed inverse problem for estimating LATs for T1 in Fig. 5. The values used in the inverse problem are 0 = 0.5 mm, 0 = 1 mm, and 0 = 2 mm, while the true value used in the forward model for generating the electrograms is 0 = 1 mm. As can be seen, using smaller values for 0 , to some extent, will still provide better results than SD . This is because the distance kernel has a larger value on its central element and using it in the inverse problem will sharpen the electrograms.

Spatial sampling resolution
The required inter-electrode distance for spatial sampling of the tissue using the electrode array, depends on the spatial bandwidth (Nyquist rate) of the electrograms. Although the true spatial bandwidth of the underlying electrogram is at the micrometer level and unknown, it is unavoidably filtered by the response of the measurement electrodes. In other words, the blurring caused by the distance kernel 0 limits the spatial resolution of the recorded electrical activity. The distance kernel (also known as point spread function (PSF) in the image processing literature [19,20]) performs as a low-pass filter and inevitably removes the high-frequency components of the recorded signal.
The so-called ''Full Width at Half Maximum (FWHM)'' of the distance kernel is commonly used as a short-hand measure of the required spatial resolution [20]. It is equal to the width between the two points on the distance kernel where the amplitude has dropped to half of its  maximum value. Fig. 9(a) shows the symmetric 2D distance kernel of size 60 × 60 mm 2 , calculated using Eq. (3) with 0 = 1 mm and with ignoring the electrode diameter. The distance kernel is normalized to have a maximum of 1. Fig. 9(b) shows a cross section of the distance kernel at the central row where the FWHM is also denoted. For this kernel the FWHM, and therefore the maximal spatial sampling distance, equals 3.4 mm. Although sampling at Nyquist rate guarantees perfect reconstruction of the sampled signal, a more conservative approach would be to record the electrograms at smaller distances than the FWHM. Therefore, the spatial resolution we use in our clinical setup equals 2 mm. Note that, a denser electrode array will only improve the signal-to-noise ratio but will not reduce the blurring caused by the distance kernel. Also note that the estimated maximal spatial sampling distance is strongly dependent on the selected electrode height 0 . Fig. 10. The logarithm of RMSE in LAT estimation using SD and SD with different values of for the five randomly generated tissue based on the T4 pattern as in Section 3.1. The logarithm of the RMSE is normalized so that the steepest descent has an amplitude of 1.

Electrode size
The distance kernel, and therefore the measured electrograms, also depend on the diameter 0 of the electrodes. In this section, we investigate the effect of the electrode's diameter 0 on the measurement model. We describe the effect of the non-zero diameter by an averaging of all per cell electrograms that the electrode's surface covers. This can be modeled by updating the distance kernel in Eq.
where function ( , ) is 1 if √ 2 + 2 < 0 ∕2, and 0 otherwise. The convolution results in an averaging effect over the electrode surface. Fig. 9(c and d) shows the distance kernel after considering the electrode diameter, calculated using Eq. (8). In Fig. 9(c) the electrode diameter is 0 = 0.45 mm. This is equal to the electrode size we use for recording our clinical data. Fig. 9(d) also shows the PSF for 0 = 1 mm. As can be seen, the electrode size does not change the morphology of the distance kernel considerably and its effect will be ignored in this study. This is because of the large value of 0 which is inherent of the tissue and cannot be improved. This result also implies that decreasing the electrode diameter to smaller values than we currently use, will not considerably improve the resolution of the electrogram recording.

Regularization parameter
As discussed in Section 2.3, the inverse problem of transmembrane current estimation is highly ill-posed and may result in many unfavorable solutions. To overcome this issue, we added a regularization term to the optimization problem in Eq. (6). Although many regularization terms can be used to stabilize the results, we opt for the sparsity of the first-order time derivative of the transmembrane current. This helps to preserve the sharp deflections of the transmembrane current (that are later annotated as LATs) while performing deconvolution. However, one can, dependent on the application, employ different regularizations such as imposing a sparsity constraint on the transmembrane current to provide sharper results, or minimizing the error between the estimated currents and a stereotype transmembrane current.
The regularization parameter controls the importance of the regularization term. Fig. 10 shows the logarithm of the RMSE in LAT estimation using SD and SD with different values of for the five randomly generated tissue based on the T4 pattern in Fig. 1. Although using larger value for makes the resulting transmembrane current sharper, as can be seen, it increases the error in LAT estimation drastically. This is because the deconvolution error is ignored. On the other hand, decreasing the value of after a certain point, causes almost no regularization and the error will stay constant.

Clinical results
In this section, we apply the proposed method to clinically recorded data. The epicardial electrograms used in this study were recorded using an 8 × 24 unipolar electrode array ( = 192) with a 2 mm inter-electrode distance and 0 = 0.45 mm electrode diameter. The electrode array is subsequently positioned visually by the surgeon on 9 mapping atrial sites using anatomical borders. The array records 5 s of sinus rhythm and 10 s of induced atrial fibrillation (IAF) signals at each site. This technique was performed in more than 400 patients of 18 years and older, with coronary and/or structural heart disease, with paroxysmal AF, electively scheduled for cardiac surgery and a ventricular ejection fraction above 40%. The acquired signals are amplified, filtered (bandwidth 0.5 to 400 Hz), sampled (1 kHz), analogue-to-digital converted (16 bits) and recorded on a hard disk. More details on the mapping approach and the electrode array specifications can be found in [21].
As with the simulated data, we use the steepest descent of the clinical electrograms (SD ) as well as the steepest descent of estimated transmembrane currents (SD ) to estimate the activation map of the electrogram array. However, we do not have access to the true LATs to evaluate and compare the performance of these two approaches. For this reason, we manually inspect the electrogram array to explain the differences in the resulting activation maps and to identify the local deflections for the clinical examples shown in this section. Although this is a quite complicated task, some spatio-temporal criteria of the electrogram array can help us. The common criteria of nonlocal activities are: (i) having no evidence of propagation [22], (ii) having time-equivalent deflections with lower amplitudes at neighboring electrodes whose amplitudes decrease with distance [23], and (iii) matching highly with the average of the neighboring electrograms [22]. While none of these criteria alone can perfectly single out the local deflection, observing some of these criteria together increases the probability of a correct annotation.
In the following subsections, we provide four categories of common cases where our proposed approach provides different results compared to the steepest descent. In each case, we use one representative example (recorded from different patients) to visualize and discuss the results. Notice that these results are based on our visual observation of limited number of patients and a more robust validation on a larger number of patients is deferred to future work. The clinically recorded electrograms for one array (8 × 24 electrodes) and Eq. (7) are used to estimate the transmembrane currents [ , , ] and subsequently the local deflection. It takes about 1.4 s to solve the deconvolution problem (with 100 iterations) for the clinical data of size = 24 * 8 without increasing the resolution, and about 6.7 s if we increase the resolution by 9 times to have smaller modeled cells (3 times increase in each axis) on a 2.9 GHz Intel i5 MacBook Pro and MATLAB R2019a.

Turning a smooth propagation into a non-smooth propagation
The most noticeable difference between the activation maps (AMs) generated from the electrograms SD and our proposed approach using the estimated transmembrane currents SD is that, in general, the AMs estimated from electrograms tend to be smoother. An example of such a case is shown as P1 in Fig. 11. As can be seen, SD introduces some heterogeneities in the area denoted by a red circle in the map, while the SD shows a homogeneous propagation. To investigate the results, we have plotted a 3 × 3 subsection of the electrograms recorded in the area of interest. The electrode locations are denoted by white dots on the AMs. The blue vertical lines denote the LAT estimated using SD and the red vertical lines denote the LAT estimated using SD . As can be seen, the recorded electrograms are fractionated, which is an indication of a zone of slow conduction in the tissue [24]. This observation matches with the AM we get using SD , but not with the AM we get using SD . This is similar to the example already shown in Fig. 5.  Fig. 11. Each row of the figure shows a different clinically recorded electrogram array denoted by P1 to P3. The first column shows LATs estimated from the electrograms SD and the second column shows the LATs estimated from the estimated transmembrane currents SD . In each case some representative electrograms belonging to the locations denoted by red dashed circles and white dots are shown on the right hand side of the figure. The blue vertical lines denote the LAT estimated using SD and the red vertical lines denote the LAT estimated using SD . The red arrow in P3 shows a new distinct wavefront in the activation map that is only visible in the activation map estimated using SD . The plot in the right-hand side of P3 shows the electrogram (4,4) and the average electrogram of its 8 neighboring electrodes.

Changing the location of the line of collision or (functional) conduction block
Another noticeable difference between the AM generated by SD and SD , is in the location of the line of collision or the line of (functional) conduction block. As can be seen in P2 in Fig. 11, a distinct late wavefront (color coded by purple) invades the mapping area from above. The AM estimated using the SD shows that this wavefront covers more parts of the mapping area compared to the AM estimated from SD . It also shows that the two distinct discontinuities in the SD map are connected. To investigate the results, we have plotted a 3 × 3 subsection of the electrograms recorded in the area of interest. The electrode locations are denoted by white dots on the AMs. The blue vertical lines denote the LAT estimated using SD and the red vertical lines denote the LAT estimated using SD . As expected, the electrograms show double potentials, which is an indicator of a collision or a functional block [24]. We can also explain the low amplitude of the second deflection by the small area its corresponding wavefront covers. On the other hand, the other wavefront (color coded in red, yellow and green) covers a bigger area and produces stronger activities. As a result, the electrodes in this area record a nonlocal activity that is much stronger than its true local activity. Moreover, on the right-hand side of the block, all electrodes have similar LAT and there is no evidence for wavefront propagation. This is an indicator of superposition of far-field activity. However the local activity seems to be very small such that it cannot even be detected using estimated transmembrane currents.

Introduction of more distinct wavefronts in the mapping area
In some cases, the proposed approach introduces new wavefronts in the mapping area. P3 in Fig. 11 shows an example where the estimated activation map using SD contains two main distinct wavefronts, while the activation map estimated using SD contains a new wavefront color coded by blue and denoted by a red arrow. We have also plotted three electrograms whose locations are denoted by white dots on the activation map. As can be seen, each electrogram belongs to one of the observed wavefronts. SD and SD provide a similar result in LAT estimation for (4, 3) and (4, 5), while they annotate two different deflections for (4,4). As can be seen, (4, 4) is composed of 3 main deflections, two of each coincide with its two neighbors, which implies these are the nonlocal activities. However, the second deflection does not match with the two wavefronts in its neighborhood. This indicates that the second deflection belongs to the local activity. This conclusion is more clear when we compare (4,4) with the average activity of its 8 neighboring electrodes shown in the right-hand side of P3 in Fig. 11. As can be seen, the difference is more distinct around the second deflection.

Providing more consistent maps over succeeding beats
One remarkable difference between an AM generated by SD and SD , is that the proposed approach provides more consistent activation maps over succeeding atrial beats. Fig. 12 shows the activation maps generated by both approaches for three succeeding atrial beats for the same location. As can be seen, both approaches show some large conduction delays in the area denoted by the red oval. However, the pattern is more consistent in the maps estimated using the proposed approach (second row). Note that P4 (in Fig. 12) and P2 (in Fig. 11) are similar. In Section 4.2 we provided some evidence to explain the difference between the activation patterns generated by SD and SD . Fig. 13 also demonstrates two example electrograms denoted by white dots in the line of block of P5. The first-order time derivative of the electrograms and the estimated transmembrane currents are also plotted. As can be seen in both cases, the amplitude of the smaller deflection increases after the deconvolution. This indicates that this small activation is a local activity and not the far-field effect of the neighboring cells.

Discussion and conclusion
In this paper, we proposed a new approach for a better estimation of local activation times for atrial mapping by reducing the spatial blurring effect that is inherent to electrogram recordings, using deconvolution. Employing sparsity based regularization and first-order time derivatives in formulating the deconvolution problem, improved the performance of transmembrane current estimation. We solved the regularized deconvolution problem efficiently using the Split Bregman algorithm. We also discussed, in summary, the required electrode array specifications including the electrode height, the electrode size and the required spatial resolution for an appropriate deconvolution. Unlike other approaches, our algorithm uses electrophysiological models to incorporates spatial information to improve LATs. We also showed, using simulated realistic tissue, that our algorithm outperforms two reference approaches: steepest descent (SD ) as the most common approach for LAT annotation, and the approach based on the spatial gradient (SG ) which also requires multi-electrode recording. Our results on simulated and clinically recorded data, show that SD is prone to make the activation maps smoother by always selecting the steepest descent and ignoring local small deflections. This arises in situations where the local activities belong to a wavefront that only invades a small part of the mapping area. The activities generated by this wavefront are thus smaller in amplitude compared to the nonlocal activities that belong to a wavefront that covers a larger area. The deconvolution, on the other hand, reduces the effect of nonlocal activities and amplifies the local activities which results in a better estimation of true local deflection. This can result in estimation of AMs that are less smooth and more irregular with probably more distinct wavefronts. It also changes the location of estimated (functional) conduction blocks or collision lines in the tissue.
While the focus of this paper is to develop a new algorithm for a better estimation of LATs in electrogram arrays, future work will include an evaluation on the effect of these changes on the further analyses performed on LATs regarding AF. More specifically, we aim to investigate the effect of the new approach on association of conduction disorders with substrates underlying AF and on the performance of existing methods in classification of AF stages in a larger group of patients. Moreover, to provide a more robust claim on the consistency of the estimated activation map over time, we aim to provide a measure of consistency/stability for a larger group of patients.

Limitations of the proposed method
In this study, we replaced each block of 3D tissue with a modeled cell in a 2D mono-layer tissue with constant cell-to-cell distance, ignoring the tissue curvature and thickness which might affect the results. Moreover, we assume a constant electrode height of 0 = 1 mm, to compute the distance kernel while it might be subject to spatial variation and quality of contact. Due to complications of simulating 3D tissue, the performance of the proposed approach was tested on two-dimensional simulated tissue. Since we do not have access to the true local activities of clinically recorded electrograms, we used manual annotation to analyze the performance of the approach.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.