Here we present a Mapping System Evaluation Framework (MSEF) to quantitatively evaluate clinical cardiac mapping systems using computational models. Our proposed framework includes the ability to evaluate mapping system performance under: 1) various electrode configurations; 2) various heart geometries; 3) various arrhythmias; and 4) the effect of noise on system performance. MSEF allows for blind testing of cardiac mapping system performance against a ‘silver standard’ in which the transmembrane potential is known throughout the entire heart. Our framework takes advantage of the fact that the methodology for quantifying the dynamic spatial patterns of transmembrane potential throughout the heart are well-established and robust, as exemplified in hundreds of experimental (e.g., optical mapping) and numerical (e.g., computational modeling) studies.

The framework includes a “library” of pre-computed simulations incorporating a range of activation patterns including paced beats, reentry, and “focal” beats replicated via pacing. For each simulation, activation times for each node in the computational mesh are computed using the maximum upstroke velocity of each action potential. The location of reentrant beats are computed via the computation of phase maps, identifying surface phase singularities, and then computing their “center of mass” from phase singularity density maps. Each entry in the library consists of: 1) transmembrane potential at every node sampled at 1 kHz; 2) activation times at every node; 3) the location of all paced beats (including simulated focal activity); and 4) the surface location and chirality of all reentrant waves.

Table 1 provides the chronological list of steps in the overall process of evaluating a generic cardiac mapping system (MS) using the MSEF. The process includes two participants: the “User” which is most likely the MS developer and; 2) the MSEF “Administrator”.

To demonstrate the feasibility of this framework, we present a specific implementation of the approach described above in this manuscript. Due to the large number of variables identified in Step 1) above, a comprehensive assessment of MSEF is beyond the scope of this study. The implementation presented here is comprised of: 1) two simulations of electrical activity in a healthy isotropic human left atria (2 seconds duration) comprised of paced (P), reentrant (R), and focal (F) beats; 2) three virtual electrode catheters: two idealized 64 basket arrays and one high-resolution 6 × 6 array; 3) well established algorithms to compute activation times and localize reentry from high resolution transmembrane patterns; and 4) simple generic mapping system algorithms. Video movies of these two simulations are provided in the Supplementary Material.

The monodomain equation governing electrical activation and propagation in excitable tissue was solved: χ C m   ∂ V ∂ t + I i o n + I s t i m − ∇ .   σ ∇ V m = 0 (1) where V_m is the transmembrane voltage, χ = 1,400 cm⁻¹ is the surface-area-to-volume ratio, and C _m = 1.0 μF cm⁻² is the capacitance per unit area. I _ion is the ionic current computed by coupling the monodomain equation with the Nygren cell model (Nygren et al., 1998) of an adult human atrial cell; I _stim is the stimulus current imposed during S1 and S2 stimulation. The conductivity was chosen to be isotropic with a value of 0.466 mS cm⁻¹ to match the conduction velocity of human atria of 55 cm/s (McDowell et al., 2015). The monodomain equation was solved using the finite element method using the Chaste software package (Mirams et al., 2013). Simulations were run on a high-resolution computational mesh of a human left atrium derived from a commercially available computer aided design (CAD) model by Zygote Cooperation. The CAD model was imported into Tetgen (Si, 2015) and an unstructured tetrahedral mesh consisting of 1.32 million nodes and 4.6 million elements with an average edge length of 252 µm was generated. The thickness was nonuniform and derived from patient specific imaging. The partial differential equations were solved using a backward Euler discretization with timesteps of 0.1 ms for both the partial differential equations and the cell model. The transmembrane voltage of each node was saved every 1 ms and simulations were run for a total of 2 s.

Two simulations were performed and videos of these are provided in the Supplementary Material. The first simulation is comprised of a single paced beat followed by 6 beats of figure-of-eight reentry, i.e., a pair of counter rotating reentrant waves (one clockwise denoted as “R+” and one counterclockwise denoted as “R-” when viewed from the endocardium), generated via an S1-S2 stimulus protocol. The paced beat was initiated at the junction of the posterior left atrium and the left inferior pulmonary vein and is referred to as “P” and the S2 was applied in the free wall of the septum of the left atria (LA). The second simulation is simulated focal activity and was constructed to allow a direct comparison with the reentrant beats. We simultaneously paced the locations corresponding to the center of mass of R+ and R-with inter beat intervals corresponding to the reentrant cycle lengths of each of the six reentrant beats; we refer to these patterns as “F+” and “F-”. Overall we simulated 13 activation patterns (AcPs) across the two simulations: one paced beat (P), six figure-of-eight patterns with clockwise (R+) and counterclockwise (R-) activation patterns and 6 pairs of focal beats (F+ and F-).

We choose two idealized generic basket electrode geometries comprised of 64 unipolar electrodes (8 electrodes spaced 2 mm apart on 8 separate splines) with a diameter 38 mm and one idealized “grid” electrode geometry comprised of 36 electrodes aligned in a 6 × 6 grid 4 mm apart. Basket catheters expand within the heart chamber into which they are placed and the distance between each electrode and the heart surface varies depending on the electrode spacing and the endocardial geometry. We initially considered the “worst-case” as all 64 electrodes residing on a 38 mm sphere “centered” in the LA; however, this case provided meaningless results which are not presented here. We consider the “best case” by finding the 64 LA sites on the mesh that minimize the distance from the endocardial surface to each electrode (see Figure 2), and then placing the electrodes 0.5 mm from the heart surface; we refer to this case as “basket good contact” (BGC). We also consider an “intermediate” case by placing the 64 locations at the midpoint of the line connecting the point on the sphere to the nearest endocardial site (see colored lines in Figure 2A); we refer to this case as “basket poor contact” (BPC). Although a sphere was used to derive the locations for BGC and BPC, the resulting electrode locations do not lie on a sphere; as such the distance between electrodes on a spline are not constant. The 8 × 8 electrode arrays for BGC and BPC are represented in 2-D arrays labeled A1 to H8 (see Figure 2B). Finally, we studied an idealized localized electrode ‘uniform array’ (referred to as “UA”) in which a 2 cm × 2 cm square was manually placed close to the endocardial surface over R+ and then each of the 36 virtual electrodes were moved such that they were 0.5 mm from the nearest endocardial site. The 6 × 6 electrode array for UA is represented in a 2-D array labeled a1 to f6 and shown in Figure 2C. Examples of electrograms from BGC and BPC with noise added are shown in Figure 2D.

Virtual electrograms (referred to as electrograms in this manuscript) were computed for all electrode locations for the three catheter configurations for the two simulations as: ∅ e x ′ , y ′ , z ′ = ∫ − σ e ∇ V m . ∇ 1 r d x   d y   d z (2) where r = x − x ′ 2 + y − y ′ 2 + z − z ′ 2 1 / 2 (3) where ∅ e is the extracellular unipolar potential (i.e., electrogram), ∇ V m is the spatial gradient of V m , σ e = 7 m/cm is the extracellular conductivity, r is the distance from a “source” point x , y , z within the heart to the electrode location, x ′ , y ′ , z ′ , and the integral is over the myocardium. This computation ignores the size of the electrode assuming it is a point. The integral was computed by summing the volume integrals over each element in the finite element mesh, calculated using Gaussian quadrature and using the finite element solution for V m (linear in each element).

The value of activation times (AcTs) for V m were computed as the time of maximum derivative of V m provided it was greater than a threshold value (α) with the constraint that two activations did not occur within a specific interval (β). The values of α and β were selected based on a sensitivity analysis performed on seven V m sites from the first simulation at five locations within the reentrant isthmus and two sites outside the isthmus. Specifically, we computed AcTs for thresholds of α = 0.1, 0.2, 0.3, and 0.5 mV/ms and for β = of 25, 50, 100, 150, 200, 250, and 500 ms. For all seven sites, the number and values of AcTs were the same for thresholds of 0.1, 0.2, and 0.3 mV/ms and window sizes of 100 and 150 m. Therefore, we choose values of α = 0.2 mV/ms and β = 100 m for the computation of AcTs from V m signals (Dube et al., 2009). Interpolation between samples was not employed so the resolution of AcTs was 1 ms.

The algorithm for identifying reentrant patterns for the simulations involved computing the 3-D filaments using state-space phase analysis using Eqn (Hosseini et al., 2017). as described previously (Pathmanathan and Gray, 2015; Galappaththige et al., 2019). θ t = atan ⁡ 2 V m t + 2 + 30 , V m t − 2 + 30 (4) where θ is a computed phase variable; endocardial phase singularities (PSs) were computed from the proper end of these filaments. PS density maps were computed using a custom Python script that calculated the number of times a PS occurred at each node within the simulation interval between 1 and 2 s. Since we identified a relatively stable figure-of-eight reentrant pattern via visual inspection, a k-means clustering algorithm was implemented to identify two clusters corresponding to clockwise (+) and counterclockwise (−) reentry on the endocardial surface. The center of mass of these two clusters were considered the locations of the two endocardial surface PS locations, R+, and R-.

The algorithm for computing AcTs from the electrograms was identical to that used for V m with the exception that the sign of the ‘derivative threshold’ was opposite; the value of ß was 100 m and the value of a was −1 mV/ms for BGC and a = −0.5 mV/ms for BPC and UA. Due to the significant differences between the morphology of V m and ∅ e signals, we did not employ phase analysis for the algorithms to identify activation patterns from the electrogram data. Instead we developed very simple algorithms to identify focal (F) and reentry (R) patterns using only AcTs. The algorithms include two parameters (a ‘difference threshold’, δ in ms; and an interval, γ in ms). We identified the presence of both F and R patterns at each electrode location using the value of AcT at that site and the AcTs of the eight surrounding electrode neighbors. A site was classified as F if all the AcT differences of the 8 neighbors and the central pixel were between -δ and γ+δ. A site was classified as R+ (R−) if there was a clockwise (counterclockwise) progression of AcT’s along the path of the 8 neighboring electrodes including a continuation of activation between beats; specifically, each of the differences along the path had to be between -δ and γ+δ. These pattern identification algorithms include the computation of eight differences and we chose δ = 2 and γ = 100.

Noise was included by adding uniformly distributed random values to ∅ e . We choose the level of noise to be equal to the threshold value (a) which varied with electrode configuration (see above) which is a level at or above clinical values (Unger et al., 2019) although the effect of noise (as a factor of threshold) is included in the Supplement. Recall that in the actual use of our proposed framework the user will provide information on the actual level of noise for their MS to the MSEF.

Evaluation of AcTs was performed for the electrograms for both BGC and BPC MSs by comparing to the corresponding values computed from V m . The acceptable level of difference in AcTs between the MSs and simulation is unclear and may depend on the activation pattern, therefore we introduce a “tolerance” variable (Tol) and analyzed performance as a function of Tol. The ability of each MS to identify AcTs was computed by identifying: 1) correctly identified AcTs; and 2) spurious AcTs. In addition, the average RMS of all correctly identified ACTs was computed. An “AcT Performance Metric” (AcTPM) was computed to assess the ability of a MS to identify AcTs: A c T P M = f C   ∗ 1 − f S   ∗   100 (5) where fC and fS are the fraction of correct and spurious AcTs, respectively. A value of 100 indicates perfect performance. Specifically, fC is computed as the number of AcTs for simulated electrogram that are within Tol of a corresponding AcT computed from V m of the nearest endocardial site divided by the total number of V m AcTs from that site; and fS is computed as the number of AcTs for a simulated electrogram that are not within Tol of a corresponding AcT computed from V m of the nearest endocardial site divided by the total number of V m AcTs from that site (if fS is >1, then fS is set equal to 1).

Activation patterns (AcPs) were computed at each site using the AcTs from the 3 × 3 array neighborhood and analyzed similarly and were considered correct if they were localized within 100 ms and if the distance to the true (x,y,z) location in the simulations was less than 1 cm. We define the temporal localization error (E_T) as the difference between the electrogram AcT at the centralized site and the corresponding stimulation time (i.e., S1, S3, S4, S5) and the spatial localization error (E_X) is the Euclidean distance between the centralized electrode location and the site identified as R+ or R-from the V m simulations (as described above). An ‘AcP Performance Metric’ (AcPPM) was computed to assess the ability of a MS to identify AcPs: A c P P M = f C   ∗ 1 − f S   ∗ 1 − f M   ∗   100 (6) where fC and fS are the fraction of correct and spurious AcPs and fM is the fraction of “misclassifications” defined as a wrong pattern type for a beat (matching the identification criteria for the above temporal and spatial distances of a different pattern). Specifically, a correct AcP from a virtual electrode array was defined as the identification of the identical pattern for the same beat for the “ground truth” (25 patterns: 1P, 6 R+, 6 R-, 6 F+, and 6 F-); and a spurious AcP from a virtual electrode array was defined as the identification of a pattern that did not correspond to the ground truth. AcTs and AcPs for BGC and BPC were compared to those computed for the simulation results for pacing, focal and reentrant patterns separately. In addition, the average temporal (E_T) and spatial (E_X) localization errors of the correctly identified and spurious AcPs were computed. Since the definition of the correct identification of a pattern depends on a 1 cm ‘tolerance’ the value of E_X is constrained (E_X ≤ 1).

Figure 3 illustrates the initiation of the paced beat (Panel A) and the location of the S2 stimulus which was applied 390 m after the paced stimuli (Panel B). A snapshot of activity from 6 views is shown in Figure 4 illustrating the figure-of-eight reentrant pattern. A video of the simulation is provided as a Supplementary Material. The figure-of-eight reentrant patterns from this simulation and focal patterns from the second simulation on the endocardial surface are shown in Figure 5 with the computed centers of mass of R+ and R-displayed as grey spheres. The location of these patterns in relation to the LA can be ascertained by viewing Figure 2A and the isochrone maps constructed from the 8 × 8 grid of electrodes for BGC for paced, reentry and focal activity in Figure 6. Similarly, isochrone maps constructed from the 6 × 6 grid of electrodes for UA are shown in Figure 7.

The average RMS value of the difference of AcTs computed from V m and ∅ e signals was less than 1 ms for all beats for all values of Tol ranging from 0 to 100 ms for both BGC and BPC. In fact, all values were below 0.24 ms except for the paced beat for BPC, for which the average RMS was between 0.38 and 0.58 ms. The fact that all values were less than 1 ms motivated the development of the novel performance metrics presented in the Methods Section. AcTPM and %S values are shown for BGC and BPC as a function of Tol for P, R, and F patterns in Figure 8. As expected AcTPM was always larger, and % S was always smaller, for BGC compared to BPC for all activation patterns. The trend was for AcTPM values to increase and % S values to decrease as Tol increased and reach plateau values with these values being less for BPC compared to BGC, and highest for P as compared to R and F for both BGC and BPC. The plateau for BGC was reached at Tol = 2 ms, where the corresponding value for BPC was Tol ≈ 10 ms. For BGC, the values for AcTPM for Tol = 2 ms were 100, 93.5 and 93.6 for P, R, and F respectively and the corresponding values for %S were 0.2, 2.5, and 1.8. For BPC, the values for AcTPM for Tol = 10 m were 84.9, 79.2 and 81.5 for P, R, and F respectively and the corresponding values for % S were 4.7, 2.0, and 3.9. Although the average RMS values were always less than 1 ms, the RMS SD was a function of Tol and was much greater for BPC compared to GC (Figure 8C).

The ability of the MSs (BGC, BPC, and UA) to identify the one paced (P) beat, the twelve reentrant (6 R+; 6 R-) patterns, and the 12 focal patterns (6 F+; 6 F-) are presented in Table 2. This comparison was carried out for two values of δ (2 and 10 ms) which corresponds to the values for which Tol reached plateau values for BGC and BPC respectively. The P beat was not identified for any MS (hence AcPPM = 0) with one misclassification and one spurious patterns evident for BGC only. Only 1 of 6 R-beats were identified for BGC (with 0 and 5 spurious patterns for δ = 2 and 10 m, respectively). Four (δ = 2 ms) or five (δ = 10 ms) of 6 R+ beats were identified by BGC and 2 of 6 for BPC (δ = 10 ms) while 5/6 were identified for UA; UA resulted in no spurious patterns while there were 0 for both BGC and BPC (δ = 2 and 10 m); the only misclassifications of R+ occurred for BGC, δ = 10 ms. Focal beats were identified with temporal error less than 10 ms for BGC, BPC, and UA, although only F+ beats were identified for UA (which is consistent with its placement, see Figure 2C). Overall, F beats were easier to identify than R beats for our simplified algorithm.

Two factors that affect the ability of a MS to identify patterns on the heart surface are: 1) the distance of the electrodes from the heart surface; and 2) the density of the surface projection of the MS electrode sites. These two values for each electrode are shown in Table 3 for BPC; the first number is the distance of the electrode to the nearest endocardial mesh node, and the second number is the average distance to the eight nearest projected endocardial sites. The fact that the location for nearest electrode for R+, F+ (E2) was closer and had a higher surface density compared to the location of the nearest electrode for R-, F- (D4) is consistent with the trend of better identification of + patterns sites compared to—patterns. To demonstrate the effect of these factors, the pattern identification algorithm described above was applied to the AcTs computed from the noiseless V m signals. These values of AcPPM for V m data were: P: 0 (1S); R+: 69; R-:14; F+: 83; F-:100 for δ = 2 m. These values are similar to BGC (δ = 2 ms) suggesting that BGC performed nearly as well as could be expected (except for F- suggesting that the optimal electrode resolution might fall between E2 and D4, see Table 3). The actual x,y,z location of R+ and F+ beats was located 0.52 cm from the nearest surface site which corresponded to B7. The fact that electrode B7 corresponds to a high surface projection density and a small electrode to surface distance (0.18 cm) explains why F+ was the only activation pattern identified by BPC (δ = 2 ms).