Conduction velocity vector field reconstruction and divergence analysis were applied to clinical electrograms, retrospectively available from one patient with persistent AF, who underwent a pre-ablation electrophysiological study. The study was approved by the local Ethical Committee and performed in accordance with the principles outlined in the Declaration of Helsinki. The patient gave written informed consent. During the electrophysiological study, a 20 pole PentaRay mapping catheter (Biosense Webster, Inc., Diamond Bar, CA, United States), composed of five radiating splines, each carrying four electrodes was sequentially moved in the LA. Twenty-one atrial regions were mapped in the patient, sampling the PVs and LA body areas. Three hundred and fifteen atrial electrograms (i.e., 15 bipolar electrograms × 21 sites) of 2 s length were recorded during the study and exported for off-line analysis. Electrograms with inadequate signal-to-noise ratio were excluded from subsequent analysis. Activation time series were automatically extracted from each bipolar electrogram as previously reported (Faes et al., 2002; Masè et al., 2015). Briefly, electrograms were pre-processed to remove ventricular interference, local atrial activation waves were identified by signal filtering and adaptive threshold crossing (Faes et al., 2002; Masè et al., 2015) and atrial activation times were estimated by measuring the barycenter of local activation waves (Faes et al., 2002).

Data are expressed as mean ± standard deviation (SD) or median [interquartile range (IQR)], as appropriate. Divergence values are given as median, maximal, minimal, and/or range values, as appropriate.

Figure 3 displays RBF reconstructions (top panels) and superimposed normalized CV fields (arrows) corresponding to the simulated patterns in Figure 1. The angular properties of the fields are quantified by divergence maps (bottom panels). All propagation patterns were precisely reconstructed by RBF interpolation. CV vectors, which were analytically determined by RBF approach, clearly indicated the direction of wavefront propagation, being orthogonal to isochronal lines. CV values estimated from RBF reconstructions approximated well the set values, resulting of 60.4, 54.2, 59.4, and 56.6 cm/s for patterns (a) to (d). Propagation pattern properties were quantitatively distinguished in terms of divergence analysis. Indeed, planar wave propagation (a) was characterized by almost-zero values of the divergence [range (−7.5⋅10^–3, 7.5⋅10^–3 mm^–1), in green]. Focal sites (b and c), acting as sources of the field, were marked by maximal positive divergence values (D_max = 10 mm^–1, in red) versus the almost-zero values of the surrounding area (D_median = 0.125 mm^–1). The collision line, acting as a field sink, displayed negative divergence values (D_min = −5.5 mm^–1, in blue, versus D_median = 6.25⋅10^–3 mm^–1).

The results of the stability analysis are summarized in Tables 1, 2, where reconstruction errors of CV vector magnitudes and directions are reported for the four patterns at the progressive removal of mapping sites (Table 1) and at increasing noise in activation time detection (Table 2). Reliable estimations of CV magnitudes and directions were obtained even with reduced electrogram sets (Table 1), although the number of sites necessary for the reconstructions increased with pattern complexity. The reconstruction of the plane wave pattern was not affected by the progressive removal of the electrograms. Focal patterns in homogeneous/heterogeneous tissues were reconstructed from a minimum of nine sites (six sites removed) with pointwise CV magnitude errors ∼12% and direction errors of ∼0.12 rad. Wavefront collision pattern were reconstructed from a minimum of 11 sites (four sites removed) with magnitude errors of ∼14% and direction errors of ∼0.13 rad. Median CV estimates were less affected than pointwise velocities by the reduction of sites, with percentage errors of 0.2 and −1.4% for focal and of 3.5% for collision patterns.

Activation time misdetection affected the reconstruction of all the patterns (Table 2), with a progressive increase of the estimation errors at increasing levels of noise. Noise had more severe effects on the estimation of CV vector magnitudes than on median CV estimates. Absolute errors on CV vector magnitudes for focal source patterns raised from 7.6 and 7.1% at 0% noise to 32.3 and 33.9% at 10% noise, and for wavefront collision patterns errors raised from 11.4 to 28.8%. Median CV values underestimated true CV values at a progressively higher extent with increasing noise levels. At 10% noise, median CV estimates decreased to 50.5 and 52.2 cm/s (error of −7.5 and −12.8%) in focal patterns and to 49.6 cm/s (error of −11.1%) in the wavefront collision pattern. In terms of CV vector directions, angle errors for the simulated patterns increased from a range of 0.002–0.09 rad at 0% noise levels to 0.27–0.43 rad at 10% noise. Noise effects on vector magnitudes and directions were significantly reduced by averaging CV values over few beats. At 10% noise amplitude, a 10-beat average reduced CV magnitude errors to 14–16% for focal patterns, and to 18% for the collision pattern, keeping direction estimation errors <0.18 rad in all patterns.

The precision of divergence analysis to locate focal sources in tissues with homogeneous and heterogeneous conduction properties is reported in Figure 4 for changing number of recording sites (left) and noise levels (right). The localization strategy was stable against a reduction in the number of recording sites. Accurate identification was maintained with a minimum of nine sites available (i.e., seven sites removed). Focal drivers were precisely localized from single-beat divergence maps in presence of mild levels of noise (<11 and <10%, for homogeneous and heterogeneous conduction properties, respectively). Accurate identification at higher levels of noise (<17 and <14%) could be accomplished by averaging divergence maps over 10 beats (gray lines).

The capability of the RBF framework to locate focal sources in a realistic, but controlled AF context, was evaluated by analyzing synthetic AF electrograms (Figures 5A, 6A, 7A), generated by a detailed LA model. Figures 5B, 6B, 7B show the activation and divergence maps obtained by moving the catheter from the region dominated by the focal driver (Figure 5) to the region with prevailing multiple-wavelet propagation (Figure 7).

In Figure 5 the presence of the focal driver at the center of the mapping system resulted in a centrifugal sequence of activation from internal to external recording points (i.e., from A₃₄ to A₁₂ and from C₁₁₁₂ to C₉₁₀). Focal activation was apparent from the reconstructed activation maps (Figure 5B, top) and was accompanied by high positive values in the divergence maps (e.g., at time T1, D_max = 3.91 mm^–1, in red, versus D_median = 0.10 mm^–1). The regularity of the focal pattern could be observed comparing successive single-beat activation and divergence maps and resulted in a consistent average divergence map (D_max = 3.45 mm^–1, in red, versus D_median = 0.10 mm^–1).

Figure 6 displays signals and maps from an intermediate region, with the focal driver located at the top right corner of the mapping system. Here, the presence of the focal driver was less apparent from visual inspection of the recorded signals, but it was revealed by single-beat and average divergence maps, displaying maximal positive values at the focal source (e.g., in the average map, D_max = 3.44 mm^–1, in red, versus D_median = 0.06 mm^–1). Single-beat activation and divergence maps showed that the area was invaded by wavefronts from the multiple wavelet region, which collided with wavefronts originating from the focal driver. The presence of collision lines resulted in negative divergence values (e.g., at time T4, D_min = −3.02 mm^–1, in blue, versus D_median = 0.07 mm^–1).

Figure 7 shows the complex propagation patterns observed in the multiple wavelet region. The activation sequence of the simulated electrograms suggested that the area was activated by wavefronts of changing directions. This was apparent in the beat-to-beat activation and divergence maps, which showed the presence of collision lines marked by minimal negative divergence values (e.g., at time T1, D_min = −4.40 mm^–1, in blue, versus D_median = −0.07 mm^–1). The irregularity of the patterns and the changing position of collision lines resulted in an average divergence map with almost-zero values [range = (−0.19, 0.08) mm^–1, in green].

The methodology was applied to multipolar catheter electrograms acquired in the LA of a patient with persistent AF. Three representative examples of activation and divergence maps observed in the patient in different LA regions are displayed in Figure 8. The observed patterns can be directly compared to the simulated maps of Figures 5–7.

Figure 8A displays the activation and divergence maps reconstructed from the mapping of the LA appendage region. The regular activation sequence from the internal to the external bipoles (i.e., from A₃₄ to A₁₂ and from C₁₁₁₂ to C₉₁₀) and the repetitive morphology of the signals suggested the presence of a stable focal activation pattern. The reconstructed map identified the focal site at the center of the mapping area. Wavefronts propagated centrifugally from the focal site with an average CV of 44.0 ± 4.3 cm/s and activated the region at a mean cycle length of 147 ms. The stability of the focal pattern was testified by the average divergence map, which displayed a maximal positive value at the center of the mapping area (D_max = 1.03 mm^–1 versus D_median = 0.08 mm^–1 in the surrounding area).

The mapping of the right superior PV ostium in Figure 8B showed a regular, but more complex pattern. Activation and divergence maps suggested the presence of a focal activation site at the top corner of the mapping area (D_max = 0.78 mm^–1), firing at a cycle length of 150 ± 16 ms. Wavefronts from the focal site activated the mapping area from top to bottom (i.e., from A₁₂ to A₃₄ to C₁₁₁₂) at an average CV of 49.8 ± 8.7 cm/s and collided with wavefronts from the bottom. Collision lines were suggested by the fragmentation of electrogram C₁₁₁₂ (Figure 8B) and were marked by negative values in the average divergence map (D_min = −0.82 mm^–1).

Maps from the LA floor (Figure 8C) evidenced a complex activation process. Single-beat activation/divergence maps suggested that the region was invaded by colliding wavefronts, propagating with an average CV of 48.0 ± 4.9 cm/s. The region was characterized by an average cycle length of 160 ms and higher variability (cycle length SD of 14 ms). The average divergence map displayed almost-zero values [range = (−0.24, 0.11) mm^–1], reflecting the instability of the activation process and the absence of a prevalent propagation pattern.

Overall, the patient’s mapping data showed an average LA cycle length of 150.1 ± 7.5 ms, with the fastest activity (136.9 ± 2.7 ms) recorded in the region of the vein of Marshall and the slowest (159.6 ± 2.9 ms) on the LA floor. Reconstructed CV vector fields showed that CV values in the LA ranged from 40.8 to 58.5 cm/s, with a mean value of 47.2 ± 4.5 cm/s. Beat-averaged divergence maps evidenced the presence of focal activation patterns in the region of the LA appendage, right superior PV and vein of Marshall, where maximal divergence values were observed (D_max = 1.00 ± 0.28 mm^–1). Collision lines were observed in proximity of focal sites and in ostial regions, where average divergence maps displayed minimal negative values (D_min = −1.03 ± 0.64 mm^–1). Differently, mapping sites on the LA body were prevalently characterized by complex and variable propagation patterns with more uniform divergence maps.

Our approach was based on a RBF reconstruction of activation patterns and corresponding CV vector fields in the mapping area. The RBF approach presents several features, which makes it suitable for integration with clinically available mapping systems. RBF interpolation does not require any assumptions on the spacing and/or density of the interpolation points (Fornefett et al., 2001; Kybic et al., 2002a,b). This allows integration with different clinically available mapping systems, with respect to other approaches (Rogers et al., 1997; Bayly et al., 1998; Alcaine et al., 2014; Zeemering et al., 2020) proposed in the experimental setting, which instead require regularly spaced and/or high-density latency data. In this study we demonstrated the capability of RBFs to accurately reconstruct CV fields from the analysis of simultaneous electrograms from multipolar catheters and we showed that the application of operators, such as the divergence, to the calculated CV fields could be used to identify focal drivers in the presence of complex propagation patterns. This extends our preliminary work (Masé and Ravelli, 2010), where we suggested the possibility of using RBFs to reconstruct activation patterns and CV fields from scattered latency data, consecutively acquired by electro-anatomic mapping system, during atrial pacing. In addition, in the present work we corroborated the stability of the reconstructions at a progressive reduction of the mapping sites, suggesting that the method may be able to cope with partial information loss due to inappropriate deployment and/or poor electrode contact with the endocardial surface.

A second advantage of RBF interpolation is that the methodology displays a certain degree of flexibility in reconstructing continuous activation patterns that can be present during atrial arrhythmias, such as focal activation, multiple wavelet propagation, and wave collision. This may represent an advantage with respect to previously proposed algorithms, such as cosine or ellipse model fitting (Weber et al., 2010, 2011; Roney et al., 2019) and triangulation approaches (Kojodjojo et al., 2006; Ravelli et al., 2011), which display accuracy in the estimation of wavefront speed, direction, and conduction anisotropy in the presence of a single propagating wavefront, but are not able to operate in the presence of other propagation patterns. Flexibility in reproducing different activation patterns, such as focal activation and wave collision, has been recently demonstrated by physically informed neural network, which may represent a promising approach also to quantify the epistemic uncertainty associated with these predictions (Sahli Costabal et al., 2020). Despite the flexibility of RBFs to reproduce different propagation patterns, it should be noticed that the definition of the interpolant function as a sum of continuous functions makes RBFs incapable to accurately reconstruct patterns where discontinuities and/or abrupt changes in activation time are present. Thus, although capable to reconstruct wavefronts with different curvature, RBFs may not be suitable to trace the head-meet-tail region and phase singularity of rotors, where discontinuities in phase values are present. This is exemplified in Figure 9, which displays the RBF reconstruction of a transient rotor observed in the complex activity region of the simulated AF. In the displayed time window, the electrical activity in the PentaRay mapped area was characterized by a rotational wave activating the tissue in a counterclockwise direction (Figure 9A). The RBF reconstruction (Figure 9B) was able to track the sequence of activation, showing a wave entering from the upper left region and turning to the right, but it could not reliably map the head-meet-tail part of the reentrant circuit. Given this limitation, specific techniques available in the literature, which detect rotors by examining the characteristics of the electrograms obtained from catheters (Roney et al., 2017; Ganesan et al., 2018, 2019; Orozco-Duque et al., 2019; Li et al., 2020), should be used to supplement our approach when the aim is the precise localization of rotors’ critical sites. Another pattern that may not be correctly reproduced by our method is the occurrence of conduction block. Indeed, our reconstruction algorithm assumes that activation times are correctly identified and aligned per beat at different atrial sites. In the presence of conduction block and missed activations at some space locations, direct RBFs interpolation may wrongly extrapolate a continuous electrical activation in these areas. To address this problem, in future implementations, missed activation times should be properly marked on the activation map and restrictions to interpolate activation patterns in these regions may be posed. Alternatively, imposition of late activation times at these sites may be considered to mimic conduction blocks in terms of extremely slow conduction areas.

Radial basis functions provide an analytical formulation of the activation patterns, which can be directly used for the analytical determination of the CV vector field and its properties in the mapping area. CV fields integrate and complete the information content of activation maps, providing quantitative data on local CVs and directions of propagating wavefronts. In combination with spatial and/or anatomical information, CV maps evidence spatial heterogeneities in conduction and slow conduction areas. As well, specific angular properties of CV vectors mark areas of focal activation and/or wave collision (Fitzgerald et al., 2005). Of note, in a study performing a quantitative comparison of the use of vector maps and isochrone cardiac activation maps to identify patterns and features associated with arrhythmias, the former displayed superior performance in mapping simple arrhythmias, reducing the number of measurements necessary to select the correct ablation target and presenting more rapid learning curves (Fitzgerald et al., 2004). Despite the potentiality of CV vector field representation of activation patterns (Fitzgerald et al., 2004), CV field maps have been mostly limited to experimental models (Bayly et al., 1998; Eijsbouts et al., 2003; Fitzgerald et al., 2003) and open-heart surgery settings (Hansson et al., 1998), and are only recently entering electrophysiological mapping data software (O’Shea et al., 2019; Williams et al., 2021). Our approach makes the construction of a CV vector maps equivalent to that of an activation map, in terms of both clinician skill and time expenditure, which supports clinical applications. Indeed, the analytical formulation derived from RBFs allows the automatic determination of the direction of propagation (given by the gradient of the function), solving the problem of precisely delineating the activation front. The interpolation function is globally estimated from all latency data, without need of point latency grouping as required by triangulation and bilinear fit model approaches (Fitzgerald et al., 2003; Kojodjojo et al., 2006; Ravelli et al., 2011). The potential of RBF approaches to measure CV is testified by a recent study, where RBFs were chosen to estimate CV values in an open software for mapping data analysis (Williams et al., 2021). Together with the computation of CV fields by RBFs, the present study implemented the application of operators, such as the divergence, to the CV fields. We showed the divergence able to reveal and quantify important arrhythmia features, such as the presence of focal sources or wave collision lines. These features may thus be used to inform the mapping strategy in order to accelerate the detection of these critical sites in the absence of human pattern recognition. Few previous studies have proposed algorithms to navigate the mapping catheter toward focal or early activation sites (Roney et al., 2014; Weber et al., 2017; Ganesan et al., 2019). An algorithm based on iterative regression analyses displayed high accuracy to predict the earliest activation site during focal tachycardias, requiring a significantly lower number of mapping points with respect to an operator-guided mapping (Weber et al., 2017). However, the method was developed to work in combination with a point-by-point acquisition mapping and thus with a single bipole at each catheter location, while instead our approach has the advantage of condensing in the divergence index the collective information from a multipolar catheter at each mapping location. Collective information from multipolar mapping systems was extracted in another study (Roney et al., 2014) by fitting a single planar or circular wave model to the activation times of the bipoles. The direction of the estimated wave was used to guide the catheter toward the earliest activation focal source. As our approach, this method had the advantage not to be limited to a specific catheter configuration, although better detection of focal sources was obtained in combination with PentaRay or spiral catheters than with circular catheters. While the approach in Roney et al. (2014) provided a global evaluation of the direction of propagation in the mapped area by a single-wave fit, our method allows a more local and precise description of the actual wave propagation patterns. In the absence of a source in the mapped area, local RBF-estimated CV vectors may be nonetheless easily combined in a macroscopic average vector to determine a prevalent propagation direction and thus direct the mapping process. More recently, a multiparameter algorithm was proposed to iteratively guide the incremental movement of circular catheter through the atrium until source (either focal or reentrant) detection (Ganesan et al., 2019). The detection was based on the computation of four indices describing the propagation pattern and the fulfillment of source detection criteria. In particular, focal sources were detected based on the “so-called” wave divergence index, which was approximated as the SD of the direction of the velocity vectors computed through subgrouping of the mapping points into triads (Ganesan et al., 2019). Differently from this method, our divergence approach may have the advantage of not requiring any subgroup of electrodes nor specific catheter configurations or symmetries, and of providing a divergence evaluation directly based on a physical definition. Our method may be conveniently integrated into these mapping navigation schemes, where the RBF-based detection of focal sources would complement the rules for the identification of rotors (Ganesan et al., 2019). However, the implementation would require the estimation of reliable threshold values for the divergence index, which may be determined by receiver-operating characteristic analysis or statistical approaches (Schreiber and Schmitz, 2000).

It is important to notice that our approach shares with all these methods the necessity of a correct detection of activation times. Indeed, although the use of an interpolation approach allows flexibility for pattern reconstruction, interpolation is more sensitive to activation time misdetections and noise effects, with respect to techniques based on model fitting (Bayly et al., 1998; Fitzgerald et al., 2003; Weber et al., 2010; Alcaine et al., 2014) or approaches that do not require activation time detection (Richter et al., 2011; Rodrigo et al., 2016; Alcaine et al., 2017; Luengo et al., 2019; Handa et al., 2020). In order to reduce inaccuracy in activation time estimations, activation waveforms from patient data were automatically identified by a well-established technique (Botteron and Smith, 1995; Faes et al., 2002), and activation times were set at the waveform barycenter (Faes et al., 2002; Masè et al., 2005, 2015). As suggested in several works (Holm et al., 1996; Faes et al., 2002; El Haddad et al., 2013; Ravelli and Masè, 2014), the use of a morphology-based activation detection, such as the barycenter method, improves estimation accuracy in the presence of fragmented electrograms. In addition, the barycenter identifies the central point of the activation waveform and thus can be spatially associated with the midpoint of the bipole, where bipolar electrogram coordinates were set. The importance of an accurate estimation of activation series was pointed out by our computer simulations, which showed a progressive deterioration of estimation accuracy at increasing noise levels. Noise affected at a higher extent the estimation of CV vector magnitudes, while vector directions (and thus divergence analysis) were less affected. On the other hand, median CV displayed a decreasing trend with increasing noise, which may be partially related to a minimization of the wavefront curvature performed by RBF. Since simulations showed the method to be more robust against electrogram loss than misdetection, it should be preferable to exclude activation series from noisy or excessively fragmented electrograms, although extrapolation of activation patterns to poorly mapped regions should be considered with caution. Alternatively, uncertainty in local activation times may be profitably addressed using novel approaches based on probabilistic interpolation, able to keep into account activation time errors at different sites and to pin the activation map more strongly at site with higher precision measurements (Coveney et al., 2020). Consistently with previous works (O’Shea et al., 2019), our simulations also suggested that, in the presence of stable patterns, noise effects on CV estimation and source localization could be reduced by averaging values over few beats. Thus, although the method is potentially able to work on a beat-to-beat basis, sequential divergence maps should be computed to distinguish transient instances of focal activation [e.g., due to epicardial breakthroughs (de Groot et al., 2010)] from the presence of a stable localized source, whose activity should be repetitive and persist over longer periods (Takahashi et al., 2006; de Groot et al., 2010; Ravelli et al., 2012, 2014). Assuming an atrial cycle length <200 ms during AF, averaging of 10 beats would require very short (<2 s) signal windows, which are consistent with clinical mapping times. The method may thus be used to complement other descriptors of the fibrillatory patterns, such as causality-based approaches (Richter et al., 2011; Rodrigo et al., 2016; Alcaine et al., 2017; Luengo et al., 2019; Handa et al., 2020), which may require longer signal windows for the analysis.

The proof-of-concept application of divergence analysis to clinical multipolar electrograms revealed the presence of stable focal activation patterns at the PVs, vein of Marshall, and LA appendage during AF. The position of the detected focal sites is consistent with previous studies in AF patients (Haissaguerre et al., 1998; Hwang et al., 2000; Schmitt et al., 2002; Nitta et al., 2004; Sanders et al., 2005; Ravelli et al., 2012, 2014). Arrhythmic episodes of multifocal origin, with triggers located outside the PV region, were observed in AF patients by multisite biatrial mapping using a basket catheter or a non-contact mapping system (Schmitt et al., 2002). Rapid repetitive activity from the LA veins, including the PVs (Haissaguerre et al., 1998) and the vein of Marshall (Hwang et al., 2000), were reported to trigger paroxysmal AF. Dominant frequency analysis applied to atrial electrograms in paroxysmal AF showed that the PVs and ostial regions were most likely to harbor a high frequency source (42%), while the probability decreased in other atrial regions and in the coronary sinus (Sanders et al., 2005). In patients with permanent AF and mitral valve disease (Nitta et al., 2004), intraoperative mapping of the entire atrial epicardium showed LA focal activation from the posterior region adjacent to the PVs and the LA appendage. Consistently, electro-anatomic mapping and combined cycle length/wave similarity analysis in patients with persistent AF localized AF sources in the PV region in 47% of patients and in the LA appendage area in 12% of patients (Ravelli et al., 2012). Our analysis revealed the presence of collision areas and more complex propagation patterns in proximity of the focal sites, as observed at the right superior PV (see Figure 8B) and in other ostial regions. This is consistent with experimental (Kalifa et al., 2006) and clinical studies (Stiles et al., 2008; Ravelli et al., 2014; Kochhäuser et al., 2017), which identified zones characterized by propagation pattern variability and fractionated activity alongside areas of fast and regular atrial activation.

In the present study the analysis of clinical data was restricted to a proof-of-concept, being limited to the LA mapping of a single AF patient. Nonetheless, the consistency of our results with previous studies suggests the potential of the approach and fosters the performance of systematical studies on larger patient groups to investigate the spatial distribution of focal activation patterns in AF and their correlation with ablation outcomes. In addition, the developed techniques may be useful to map more organized forms of AF, as well as secondary atrial tachycardias at re-do procedures (Haissaguerre et al., 2006). Further validation of the method in larger clinical datasets, including different types of atrial arrhythmias, is necessary to clinically validate the method and to assess its applicability and benefit for the optimization of ablation strategies.