Novel algorithm for accelerated electroanatomic mapping and prediction of earliest activation of focal cardiac arrhythmias using mathematical optimization

Background Premature beats (PBs) are a common finding in patients suffering from structural heart disease, but they can also be present in healthy individuals. Catheter ablation represents a suitable therapeutic approach. However, the exact localization of the origin can be challenging, especially in cases of low PB burden during the procedure. Objective The aim of this study was to develop an automated mapping algorithm on the basis of the hypothesis that mathematical optimization would significantly accelerate the localization of earliest activation. Methods The algorithm is based on iterative regression analyses. When acquiring local activation times (LATs) within a 3-dimensional anatomic map of the corresponding heart chamber, this algorithm is able to identify that exact position where a next LAT measurement adds maximum information about the predicted site of origin. Furthermore, on the basis of the acquired LAT measurements, the algorithm is able to predict earliest activation with high accuracy. Results A systematic retrospective analysis of the mapping performance comparing the operator with simulated search processes by the algorithm within 17 electroanatomic maps of focal spreading arrhythmias revealed a highly significant reduction of necessary LAT measurements from 55 ± 8.8 to 10 ± 0.51 (n = 17; P < .0001). Conclusion On the basis of mathematical optimization, we developed an algorithm that is able to reduce the number of LAT measurements necessary to locate the site of earliest activation. This algorithm might significantly accelerate the mapping procedure by guiding the operator to the optimal position for the next LAT measurement. Furthermore, the algorithm would be able to predict the site of origin with high accuracy early during the mapping procedure.


Introduction
Premature beats (PBs) are a common finding in patients with structural heart disease, but they can also occur in otherwise healthy individuals. In patients with drug refractory symptomatic PBs or frequent monomorphic ventricular PBs in patients with reduced left ventricular ejection fraction, catheter ablation is well indicated. 1,2 Since the introduction by Gepstein et al, 3 3-dimensional (3D) electroanatomic mapping systems are increasingly applied to locate the exact site of origin of PBs.
However, infrequent occurrence of PBs during the procedure can hinder the creation of a detailed activation map within an acceptable period of time, thereby limiting procedural success. In these particular cases, it may be reasonable to first reconstruct the anatomy of the heart chamber during sinus or paced rhythm, adding the local activation times (LATs) to the existing anatomical information within a second step as a so-called remap. Figure 1A exemplarily displays the electrocardiogram of a young patient with long QT syndrome suffering from recurrent episodes of torsade de pointes tachycardia triggered by short-coupled left ventricular PBs. In this particular case, the operator decided to first reconstruct left ventricular anatomy. By obtaining LAT measurements (exemplarily displayed in Figure 1B  Motivated by established systematic search routines as, for example, applied for the rescue of avalanche victims, we strived for developing an algorithm for optimized data acquisition to accelerate the mapping procedure in cases of rare arrhythmia occurrence. 4 Our strategy was based on the assumption that when deciding about the localization of each next LAT measurement, the amount of additional information may largely differ depending on the exact location of the measuring point. We therefore developed a mapping algorithm that is able to calculate the amount of additive value at each nodal point of the geometry and automatically position the next LAT measurement at the site of maximum additive information. Furthermore, the algorithm is able to predict earliest activation by extrapolation on the basis of the acquired LAT measurements with high accuracy.

Electrophysiological procedures and data acquisition
Seventeen patients who underwent catheter ablation of focal arrhythmias guided by a 3D mapping system (CARTO 3, Biosense Webster Inc., Diamond Bar, CA) between March 1, 2014 and August 31, 2015 were selected retrospectively from our database. The study was approved by the local ethics committee and was performed in accordance with the Declaration of Helsinki (64th WMA General Assembly, Fortaleza, Brazil, 2013). Data were recorded and analyzed using the LABSYSTEM PRO electrophysiological recording system (Boston Scientific, Marborough, MA,). Electroana-tomic maps were established with the point-by-point acquisition mode of the CARTO system.

Development of patient-specific geometric models
For each patient, the raw data of the geometry of the heart chamber exported from the CARTO system consist of a triangular mesh. Some points (measurement points) are further labeled with a LAT obtained by the operator. On the basis of the spatiotemporal information exported from the CARTO system, we established a 3D geometry for each patient using the MATLAB software package (MathWorks, Natick, MA). Figure 2A schematically displays the structure of the surface of this mesh/graph. Within this graph, the distance between 2 nodal points was measured using the shortest distance over the connecting edges. The shortest distance between the 2 nodal points i and j is denoted as d ij .
Considering the focal character of the arrhythmia, we assumed that electric excitation would spread centrifugally over the connecting edges of the mesh geometry with a constant conduction velocity (CV). CV for each patient was calculated on the basis of the acquired LATs and the geometric positions of the mapping points (Table 1). When fitting a straight line through these data, CV could be simply obtained from the slope of the line (Supplemental Figure 1). The LAT at a certain nodal point (t i ) is given by a linear model using the shortest path distances and the following formula: where t i is the arrival time of the signal at point i, d ik is the shortest distance between nodal point i and source k, CV is the conduction velocity, t 0 is the unknown earliest activation time, and SEM is an individual measurement error. This error had to be included because the estimated CV obtained from the exported nodal points exhibited a certain degree of uncertainty, reflected by the r 2 value ( Table 1).

Prediction of earliest activation
When the algorithm obtained a measurement within the previously generated geometric model, the LAT at this specific nodal point was computed by the distance to the origin and the CV using formula 1. On the basis of the location and LAT of the obtained measurements, the algorithm predicted the origin of the signal by solving a linear regression problem for every nodal point of geometry j. To estimate the CV and the initial time t 0 , we minimized the objective Solving this problem for every nodal point on the geometry, the nodal point exhibiting the best fit (highest r 2 value) was considered the origin of the signal or at least the best estimate of the origin on the basis of the available information. An example is considering that a number of mapping points  have been acquired and the algorithm tries to identify the site of earliest activation. In this case, it would go through all other nodal points each time establishing a regression analysis (see Supplemental Figure 1) with all available mapping points. Figure 2 explains the search routine within a simple planar mesh geometry. The red nodal point marks the true origin of excitation, and the black points (labeled a, b, and c) mark the measurement points already obtained at this time point. When assuming the blue point as the possible source, the calculated distances to the measurements points are 1, 3, and 3 numbers of edges. However, as the LATs at the measurement points in fact characterize the distance to the true origin (red point), there is no good correlation between distance and time ( Figure 2C). In contrast, when assuming the red point as the possible source, the distances as well as the time delays to the measurement points (black points) are 1, 2, and 3. Figure 2D displays the excellent correlation of these points, thereby identifying the red point as true origin.
Optimizing the localization of the next measurement point Using the aforementioned strategy, the algorithm is able to predict the most likely localization of the origin. However, a main goal of our work was to develop an algorithm that considerably reduces the number of measurement points. For this purpose, the algorithm needs to select that next measurement point that increases the quality of our estimate best. To identify the optimal next measurement point, the algorithm tried to minimize the standard error of the y-intercept of the regression curves for all nodal point using the following formula: where a l is the positive weight used to define how important a nodal point and its detection as source is. The second regression parameter of the regression at nodal point i is t 0i . In other words, we wanted to minimize the variance of the y-intercept for all regressions we perform by picking the next measurement point. The reason for this approach is that every point of the regression curve possesses a certain degree of uncertainty. This is reflected by the variance (see red lines in Supplemental Figure 1). In addition to the correlation coefficient, the variance of the y-intercept (at distance 0) provides information about the quality of fit. Therefore, to identify the specific nodal point where a next mapping point would best improve the quality of the map, the algorithm again went through all nodal points checking how much a measurement at this specific point would reduce the variance of the y-intercept of the regression curves at all other points. This specific nodal point that results in the maximum reduction of the variance of the y-intercept of all other nodal points was identified as the next mapping point. The mean overall computing time for the calculation of the optimal position of the next mapping point as well as the revised prediction of earliest activation was 124 ms.

Detecting the source
In order to locate the site of earliest activation, the algorithm, at a certain point in time, stops taking mapping points based on the degree of maximum additive information and starts localizing the site of origin by a direct search close to the predicted site of origin. This change of the search strategy seems reasonable, as points collected on the basis of maximum additive information are generally remote and not close to the predicted origin. To identify the point in time, when this second phase of the search routine needed to be entered, the algorithm continuously compared the quality of fit (derived from regression analysis) between the predicted site of origin and all points within 1 cm around it. As soon as the difference of the quality of fit to the surrounding points reached a certain threshold (surrounding points significantly worse), the search routine was changed to a more or less heuristic search around the predicted site of origin.

Mapping of earliest activation by the operator
Electroanatomic maps from 17 patients who underwent ablation of focal arrhythmias were selected retrospectively from our database. On average, a number of 55.1 6 8.8 (n 5 17) LAT measurements had been acquired by the operator before ablation. The geometry of the aforementioned patient with long QT syndrome is displayed exemplarily in Figure 3A. The LAT measurements obtained by the operator are displayed as color-coded points, with red representing early activation times and purple late activation times (see timescale). The projection of the search path followed by the operator is displayed. Figure 3B displays the distance between the measuring points and the origin in relationship to the corresponding LATs. In this particular case, a total number of 27 LAT measurements were obtained. The corresponding linear regression exhibited a good correlation, yielding an r 2 value of 0.86 (slope 0.55 6 0.044 ms/ mm; y-intercept 271.1 6 1.1 ms; n 5 27). Table 1 gives an overview of the correlation coefficients, the calculated CVs, and the mapped heart chambers of all patients.

Mapping of earliest activation by the algorithm
Next, we analyzed mapping performance of the developed algorithm. For this purpose, realistic models of excitation propagation had to be first established for each patient. We therefore computed the LAT for each nodal point of the mesh geometry using the electroanatomic map created by the operator, the identified origin, the calculated CV, and the individual measurement error. These geometries served as a test environment for the algorithm. At the beginning of the automated mapping procedure, the algorithm had information only about the anatomical situation derived from the mesh geometry, such as the exact localization of the nodal points and the connecting edges, which means that the algorithm was blinded to all LATs of the geometry. At the very beginning of the automated mapping procedure, the first 3 LAT measurements were chosen en bloc because the selection of only 2 points would result in a perfect linear regression fit, rendering all points equally likely to be the origin. These 3 measurements are automatically performed by the algorithm purely on the basis of the shape of the geometry. Now that the LAT and localization of the initial 3 measurement points had been obtained, the algorithm went through all nodal points of the geometry, each time calculating the distance to the 3 initial measurement points and the linear regression for distance and LAT. The nodal point with the best fit was then assumed to be the origin (predicted origin). Figure 4A exemplarily displays the situation for the previously described patient (Figure 1) at this exact point in time. The initial 3 measurement points are displayed as black points. The quality of correlation (r 2 ) of each nodal point is color coded, with blue representing strong correlation and white weak correlation. The nodal point with the best correlation (predicted origin) is highlighted as a large blue point. The red point indicates the localization of the true origin. Since only 3 LAT measurements had been performed at this point in time, the predicted origin is still remarkably displaced from the true origin. The correlation between the distance to the measurement points and the LATs for the predicted origin (large blue point) is displayed in Figure 4B.
Next, the algorithm aimed at identifying that specific nodal point that would add maximum information about the localization of the true origin. Considering that the y-intercept of the regression curve represents the predicted origin, the algorithm searched for that specific nodal point at which an LAT measurement would reduce the variance of the y-intercept for the regression curves of all nodal points most. Figure 4C displays a mesh geometry showing the additive value at each nodal point, with black representing maximum additive information and white minimum additive information. The location of the maximum additive value is marked by an arrow (Figure 4C). Two itera-tions later and now based on the information of 5 LAT measurements, the predicted origin evidently migrates toward the true origin ( Figures 4D-4F). Figure 5A displays the situation after the seventh iteration when the algorithm located the true origin by placing the last measurement in this position. Comparable to Figure 3A, the acquired measurement points are color coded, with red representing early activation times and purple late activation times. Again, the search path is displayed on the frontal and sagittal planes. In contrast to the operator (27 measurements), the algorithm identified the true origin with 7 LAT measurements. Figure 5B displays the correlation between the distance to the origin and the LAT. To allow a direct comparison between the map automatically generated by our algorithm with the original CARTO map, we created a color-coded CARTOlike activation map ( Figure 5C). Identically to the remap displayed in Figure 1B III , our map is based on the spatiotemporal information of 7 LAT measurements. A direct comparison of both maps visualizes the obvious difference in accuracy.

Comparison between the operator and the algorithm
Next, we analyzed the diagnostic performance of the algorithm within all previously generated test geometries. Figure 6A displays the mean number of iterations that were necessary to identify the origin of the arrhythmia in each patient. Overall, the mean number of LAT measurements that were needed by the algorithm to identify the origin was 10 6 0.51 (n 5 17). Compared to the algorithm, the operator, on average, took 5 times as many LAT measurements (55 6 8.8; n 5 17; P , .0001).
However, when performing a head-to-head comparison between the mapping performance of the operator and the algorithm, it has to be taken into account that the algorithm always started mapping within an existing anatomy whereas the operator in some cases had reconstructed the anatomy simultaneously while mapping activation times. For this purpose, Figure 6B compares only those particular cases in Figure 3 Evaluation of the mapping approach of the operator. A: Reconstructed mesh geometry obtained from the ablation procedure in the aforementioned patient with long QT syndrome. LAT measurement points obtained by the operator are displayed as points with color-coded activation times, with red representing early activation times and purple late activation times (see color scale). The search path from the LAT measurement to the LAT measurement followed by the operator is displayed as a projection on the frontal and sagittal planes. B: Distance between measurement points and the corresponding LATs shows a good correlation when assuming the site of best fit as the true origin (red point). LAT 5 local activation time.
which the operator was able to map activation times within a preexisting anatomy (remap). In these 10 cases, the site of origin could be identified within 11 6 0.89 LAT measurement points by the algorithm as compared to 42 6 7.0 LAT measurement points by the operator (n 5 10; P , .001).

Clinical implications
The algorithm outperformed the operator in almost every case in terms of the number of mapping points necessary to locate the site of origin, thereby pointing to shorter procedure times. In our opinion, there are 2 main reasons for this observed effect: (1) When directly comparing the number needed based on iterative linear regression analyses, our algorithm is able to calculate the degree of redundancy from each nodal point of the geometry, thereby identifying the nodal point with maximum additive information. This optimized search routine guarantees that the site of earliest activation can be located with a low number of LAT measurements. (2) Except for optimized data acquisition, our algorithm exhibits a second feature that might contribute to the significant reduction of mapping points. The linear regression analyses calculated for each nodal point are used not only to calculate the amount of additive information but also to predict the site of earliest activation. This means that LATs of any nodal point is extrapolated from the acquired LATs. In contrast, activation maps generated by the CARTO system are based on interpolation between the acquired mapping points. This difference might account for the higher quality of the activation map generated by our algorithm. As an example, when comparing the CARTO map based on the spatiotemporal information of 7 LAT measurements displayed in Figure 1B III with our map based on the same number of LAT measurements ( Figure 5C), a clear difference in the accuracy is visible. Interestingly, a quite similar approach using regression analyses has been published recently for the identification of the source of network-driven contagion phenomena such as the 2009 H1N1 influenza pandemic or the 2003 severe acute respiratory syndrome epidemic. 5 In their work the authors report that on the basis of geographical locations, arrival times of infection, and predicted traveling times, the source of the infection can be located using correlation analyses. Compared to our simple 3D electroanatomic model, these simulations are quite complex because of the inhomogeneous network structure and limited knowledge about the exact traveling times. 6 However, motivated by the good results of these studies, it might be a promising approach to further increase the degree of complexity of our models, for example, by including anatomical structures (ie, mitral or tricuspid annulus) or by including areas of scared tissue with reduced CV, for example, to allow fast activation mapping in scar-related ventricular tachycardia.

Alternative mapping techniques
Mapping of focally spreading cardiac arrhythmias is a relevant clinical topic, and a variety of strategies have been developed to accelerate localization, thereby increasing success rates. Hocini and coworkers 1,2,7 recently published promising data from a multicenter study applying a novel technique of high-resolution noninvasive mapping for the ablation of focal arrhythmia. However, it has to be taken into account that this technique requires a thoracic computed tomography scan, an electrode vest, and an additional mapping system (ECVue, CardioInsight Technologies, Inc., Cleveland, OH). A similar approach of inverse potential mapping based on the combination of a magnetic resonance  imaging scan and a body surface potential map has recently been published by Bhagirath and coworkers. 3,8 However, noninvasive mapping techniques have not yet entered clinical routine. Similar to noninvasive mapping techniques, the use of multielectrode catheters has also been proposed to accelerate endocardial mapping of focal arrhythmias. Using noncontact 4,9 or contact 1,10,11 acquisition of endocardial electrograms, these systems allow simultaneous assessment of a larger number of LATs using basket or balloon catheters. Compared to all technologies described above, our algorithm offers the advantage that it may easily be implemented into standard 3D electroanatomic mapping systems without any need for additional hardware components or preprocedural imaging modalities.

Study limitations
Because of obvious technical reasons, the diagnostic performance of the algorithm could be assessed only retrospectively using electroanatomic maps that had been exported from the CARTO system. Furthermore, operators with different degrees of expertise conducted the mapping procedures. Therefore, it can be suspected that fewer LAT measurements might have been necessary when the operator with most experience would have performed all procedures. However, in this way, the results of the study might even be more representative for a real-world situation.

Conclusion
We developed an automated mapping algorithm for the identification of the site of earliest activation within 3D electroanatomic maps. We further show that when compared to an operator, the algorithm is able to locate the site of earliest activation with a significantly lower number of LAT measurements. When integrated into an electroanatomic mapping system, this algorithm might significantly accelerate the procedure by guiding the operator to the optimal position for the next LAT measurement, thereby reducing the number of points with a high degree of redundant information.
Furthermore, the algorithm would be able to predict the site of origin with high accuracy early during the mapping procedure.