Patients referred for endocardial and epicardial EAM and ablation were enrolled. Epicardial mapping was indicated because of either recurrent ventricular tachycardia with a suspected epicardial substrate or symptomatic premature ventricular complexes with a prior failed endocardial ablation. Anti-arrhythmic drugs, except amiodarone, were discontinued for a minimum of three half-lives prior to the ablation procedure. Amiodarone was continued because of its long half-life. The study protocol was approved by the local institutional review board (University Medical Center Utrecht, Utrecht, Netherlands; protocol nr.17/628). The study was conducted according to the declaration of Helsinki and all patients gave informed consent prior to non-invasive and invasive mapping.

The workflow of the study is depicted in Figure 1. Patients underwent 67-electrode body surface potential mapping (sampling frequency 2048 Hz, Biosemi, Amsterdam, Netherlands) prior to the invasive mapping procedure and the electrode positions were captured using a 3-dimensional camera (Intel Realsense D435, Santa Clara, CA, United States) (Hoekema et al., 1999). Per patient, cardiac computed tomography (CT, Philips Healthcare, Best, Netherlands) was performed to manually create patient specific anatomical models of the ventricles with both epicardial and endocardial surfaces, ventricular blood pool, lungs and thorax. The ventricular anatomical models were supplemented with patient specific endocardial structures associated with early ventricular activation through the His-Purkinje system (e.g., the left ventricular papillary muscles and right ventricular moderator band) (Durrer et al., 1970). Electrode positions were reconstructed by registering 3-dimensional images to the thorax model. The volume conductor model was computed using the boundary element method. Conductivity values of 0.2 S/m for the thorax and ventricular tissue, 0.04 S/m for the lungs and 0.6 S/m for the blood cavities were used (Supplementary Methods).

Baseline drift and 50 Hz noise were removed from the body surface potential map signals. Per patient, five subsequent sinus rhythm complexes were selected to be analyzed in the iECG procedure. Premature ventricular complexes and sinus rhythm complexes prior to premature ventricular complexes were excluded from analysis. The root mean square of all recorded signals was used to annotate QRS onset, J-point and T-wave end. One lead from the standard 12-lead ECG was used as timing reference to allow comparison of absolute timings between iECG estimations and invasive EAM timings.

The novel iECG method has been described in more detail in the Supplementary Methods (Greensite et al., 1990; Huiskamp and Greensite, 1997; van Dam et al., 2009; van Oosterom, 2014). In short: the iECG method simulates body surface potential maps using the patient specific EDL cardiac source model, the patient specific volume conductor model and the estimated ventricular activation sequence. Nine regions containing potential foci were localized: four at the left ventricular septum, two at the base of both the posterior and anterior papillary muscles of the left ventricle, two at the right ventricular septum and one at the insertion of the moderator band at the right ventricular free wall free wall (Durrer et al., 1970). The fastest route algorithm was used to compute activation sequences emerging from these locations and combinations of foci (van Dam et al., 2009). All possible combinations of foci were tested as the initial estimation (Supplementary Methods). The activation sequence from the initial estimation with the highest correlation between the simulated body surface potential map and the recorded body surface potential map, was selected as input for the optimization step (Supplementary Methods) (van Dam et al., 2009). The optimized activation sequence was used to assign LAT to each node in the patient specific ventricular anatomical model (Figure 1 and Supplementary Figure 1).

Invasive EAM was performed under general anesthesia during sinus rhythm or atrial pacing. Ventricular paced complexes and premature ventricular complexes were excluded from analysis. Epicardial access was obtained by percutaneous subxiphoid approach (Sosa et al., 1996) and endocardial access was obtained through the right femoral vein. Access to the left ventricle was gained through a transseptal puncture, using a steerable sheath (Mobicath, Biosense-Webster Inc. Irvine, CA, United States). Anatomical coordinates, LAT maps and voltage maps were automatically created with EAM systems (Carto-3, Biosense-Webster Inc. Irvine, CA, United States or EnSite Precision, Abbott, Chicago, IL, United States) without prior integration of cardiac CT images. Endocardial and epicardial EAM was performed with multi-electrode catheters (PENTARAY^® catheter, Biosense-Webster Inc. Irvine, United States or ADVISOR^TM HD Grid mapping Catheter, Abbott, Chicago, Il, United States). Unipolar and bipolar electrograms were simultaneously recorded with standard 12-lead ECG (band pass filters 30–500 Hz, sampling frequency 1000 Hz), and one of these leads was used as timing reference for electrograms. Post-procedure, bipolar and corresponding unipolar electrograms were manually reviewed by investigators who were blinded to the information from the corresponding iECG map. LAT was set at the maximal amplitude of the bipolar signal, corresponding to maximum downslope (dV/dt) in unipolar electrograms (see Figure 2 for examples) (Cantwell et al., 2015). In case of doubt, recordings from neighboring electrograms were taken into consideration to determine LAT. Epicardial and endocardial myocardium with abnormal voltage electrograms was defined as bipolar voltage amplitude < 1.5 mV.

Anatomical coordinates with corresponding annotated LAT and bipolar voltage, obtained during EAM, were exported (MATLAB-2017a, The Mathworks Inc, Natick, United States). These anatomical coordinates were semi-automatically aligned to the CT-based ventricular anatomical model, according to anatomical landmarks (right ventricular outflow tract and the apex of the ventricles, Figure 1). Endocardial alignment was optimized using a closest point matching algorithm (Bergquist et al., 2019). Surface Laplacian interpolation was used for areas with incomplete EAM, within a distance of 10 mm. To reduce misalignment errors, invasively collected datapoints for myocardial surfaces were projected onto the nearest node of the CT-based model and all projections per node were averaged. iECG-LAT maps were referenced to the same timing reference used during the EAM procedure. Pearson’s inter-map correlation coefficient and inter-map absolute difference in milliseconds (ms) were determined for the epicardium, right ventricular endocardium and left ventricular endocardium. Breakthrough of activation was defined as nodes with the lowest LAT value, and sites of latest activation were defined as the node with the highest LAT value. Euclidian distances between sites of earliest and latest activation were determined in millimeters (mm). Myocardial conduction velocity over surfaces was calculated as the minimum positive velocity between nodes, velocities more than 3 mm/ms were excluded. A relatively high cut-off of 3 mm/ms was used to account for velocities observed in regions with a high density of Purkinje-myocardial junctions as the conduction velocity of Purkinje fibers ranges between 2 and 3 mm/ms. This cut-off was used to take into account that the electrical pulse may spread via the Purkinje fibers to the neighboring myocardial tissue instead of via the myocardial tissue itself. Ventricular activation sequences were presented in right anterior oblique, left anterior oblique and inferior views (Cosio et al., 1999).

Data were presented as mean ± standard deviation or median [interquartile range], supplemented with 95% confidence interval (CI). Continuous data were compared using (un)paired Student’s t-test or Mann–Whitney U test as appropriate. Differences between iECG-LAT maps and EAM-LAT maps were presented as absolute difference in ms for timings or absolute difference in mm for differences in sites of breakthrough, earliest activation or latest activation. iECG-LAT and EAM-LAT maps were compared with Pearson’s linear correlation and presented as correlation coefficient. Agreement between iECG and EAM-LAT timings was quantitatively compared by Bland-Altman plots. A 2-sided P-value of <0.05 was considered significant. Statistical analysis was performed in MATLAB (MATLAB 2017a, The Mathworks Inc, Natick, MA, United States).

Thirteen patients (77% males, age 48 ± 20 years) referred for epicardial and endocardial mapping and ablation of ventricular tachycardia (n = 10) or symptomatic premature ventricular complexes (n = 3) were included. Patients were diagnosed with arrhythmogenic cardiomyopathy (n = 5), dilated cardiomyopathy (n = 2), symptomatic premature ventricular complexes (n = 3), or ventricular arrhythmias after healed myocarditis (n = 3). Patients had either sinus rhythm (n = 10) or atrial pacing by an implanted permanent pacemaker (n = 3) during body surface potential recording and the EAM procedure, see Supplementary Table 1 for a summary of the included population and Supplementary Table 2 for a detailed description per included patient.

The patient cardiac anatomical models had an inter-node spatial resolution of 8 ± 1 mm. The QRS complex morphology of the recorded body surface potential maps correlated with the QRS complex morphology of the simulated body surface potential maps in the iECG procedure (correlation coefficient = 0.97 ± 0.02). The QRS morphology of the timing reference lead during EAM correlated with the timing reference lead of the recorded body surface potential map (correlation coefficient = 0.94 ± 0.02).

Epicardial EAM was performed in all patients, right ventricular endocardial EAM in 10 patients and left ventricular endocardial EAM in four patients. EAM consisted of median 4611 [3369–5633] epicardial electrograms, 910 [280–1638] right ventricular endocardial electrograms and 605 [247–1412] left ventricular endocardial electrograms. The number of annotations per square mm was 20 ± 11 for the epicardium, 10 ± 5 for the right ventricular endocardium and 8 ± 4 for the left ventricular endocardium. The percentage of EAM per surface was on average 67 [range: 48–82]% of anatomical nodes for the epicardium, 45 [range: 15–79]% for the right ventricular endocardium and 48 [range: 22–71]% for the left ventricular endocardium. The anatomical EAM model was limited to the locations where the catheter had been positioned during the EAM procedure.

Figure 2A shows an example of the comparison of iECG and EAM for LAT maps, and the comparison between earliest and latest activated nodes for both the epicardium and endocardium. The ranges between earliest and latest ventricular activation were not significantly different between iECG-LAT maps and EAM-LAT maps (111 ± 23 vs. 124 ± 39 ms, p = 0.311). The ranges of earliest and latest activation per patient are included in Supplementary Table 3. Figure 2B shows an example of the iECG and the EAM approach in a patient with a healed myocarditis with right bundle branch block. The fast and His-Purkinje mediated activation of the left ventricular myocardium is shown in contrast to the relatively slower activation of the right ventricle due to the right bundle branch block. Figure 2C shows an example of the activation pattern and epicardial electrogram annotation in a patient with arrhythmogenic cardiomyopathy. Furthermore, all LAT and voltage maps of each included patient are available as Supplementary Figure 1. The mean correlation coefficient between iECG-LAT maps and EAM-LAT maps was 0.54 ± 0.19; [95% CI:  0.49–0.59] for epicardial maps, 0.50 ± 0.27; [95% CI: 0.41–0.58] for endocardial right ventricular maps and 0.44 ± 0.29; [95% CI: 0.32–0.56] for endocardial left ventricular maps (Table 1). The moderate agreement of LAT between iECG and EAM maps is shown in Figure 3A for all included electrograms on the epicardium (R = 0.632, p < 0.001), right ventricular endocardium (R = 0.597, p < 0.001) and left ventricular endocardium (R = 0.546, p < 0.001). Figure 3B shows that a prolonged QRS duration of the included complexes did not affect correlation coefficient or absolute difference. Figure 3C suggest that a higher density of mapped electrograms per mm2 reduces the scatter of correlation coefficients. The absolute difference for epicardial LAT maps was 17.4 ± 7.2 ms; [95% CI: 15.6–19.2], for endocardial right ventricular maps 19.5 ± 7.7 ms; [95% CI: 17.2–21.7], and for endocardial left ventricular maps 27.9 ± 8.7 ms; [95% CI: 24.2–31.5]. The relation between the percentage of mapped anatomical points during EAM and the agreement for LAT values is shown in Figure 3D. The correlation coefficient between iECG-LAT maps and EAM-LAT maps was not significantly affected by the absolute number of EAM electrograms (p = 0.324), the number of electrograms with abnormal voltage (p = 0.306) or the QRS duration (p = 0.485) (see Supplementary Figure 2). However, the annotation density and the percentage of mapped anatomical points per map affected the agreement between iECG and EAM. In maps with a low annotation density or lower percentage of mapped anatomical points the correlation coefficients were low (Figures 3C,D). This may have negatively affected the observed correlation coefficients in this study because endocardial EAM was often limited to either the right ventricular or left ventricular surface. The iECG estimations were based on five QRS complexes selected from the body surface potential maps, but a Bland-Altman analysis did not result in divergent results per included QRS complex. These scatter plots and Bland-Altman plots for each included patient are available in Supplementary Figure 3.

The number of endocardial breakthrough points was similar when comparing iECG-LAT maps and EAM-LAT maps: 3.3 ± 0.8 vs. 3.2 ± 0.9 for epicardial maps, 1.8 ± 0.6 vs. 2.1 ± 0.6 for right ventricular endocardial maps and 1.8 ± 0.5 vs. 1.8 ± 1.0 for left ventricular endocardial maps (Table 1). These findings were in line with the observations of Durrer et al. and the assumptions of the iECG initial estimation (Durrer et al., 1970). Epicardial breakthrough of activation had an absolute difference between iECG-LAT maps and EAM-LAT maps of 42.1 ± 18.6 mm; [95% CI: 36.7–47.5]. For endocardial breakthrough of activation, the absolute difference was 29.9 ± 16.0 mm; [95% CI:  25.1–34.8] for the right ventricular endocardium and 31.0 ± 16.8 mm; [95% CI: 23.8–38.1] for the left ventricular endocardium. The latest activated nodes had an absolute difference between iECG and EAM of 54.1 ± 26.9 mm; [95% CI: 47.5–60.7] for epicardial maps, 46.7 ± 28.8 mm; [95% CI: 38.8–54.7] for right ventricular endocardial maps and 32.7 ± 17.2mm; [95% CI: 25.1–40.4] for left ventricular endocardial maps (Table 1). The timing of the latest activated nodes differed 12.8 [2.9–29.3] ms; [95% CI: 6.4–31.7] for epicardial maps, 15.2 [10.1–28.7] ms; [95% CI:  8.5–32.5] for right ventricular endocardial maps and 20.8 [10.4–57.3] ms; [95% CI: 12.5–71.4] for left ventricular endocardial maps. The myocardial conduction velocity was not significantly different between iECG and EAM maps for, respectively, the epicardium (1.26 ± 0.16 vs. 1.26 ± 0.20 m/s, p > 0.999), right ventricular endocardium (1.13 ± 0.09 vs.. 0.94 ± 0.17 m/s, p = 0.069) or left ventricular endocardium (1.03 ± 0.11 vs. 0.92 ± 0.07 m/s, p = 0.968).

Previous quantitative EPD-based validation studies showed higher agreement between ventricular paced complexes and EAM, compared to sinus rhythm complexes (Duchateau et al., 2019; Graham et al., 2019). Duchateau et al. (2019) showed poor epicardial inter-map correlation coefficient (−0.04 ± 0.3) during sinus rhythm, although correlation coefficients increased with increasing QRS duration. This relation is most likely explained by the complexity of multiple simultaneous ventricular activation waveforms occurring during sinus rhythm, which decreases in rhythms with a single focus (Duchateau et al., 2019). In the present study, a considerably higher agreement (correlation coefficient 0.54 ± 0.19) between EAM and the novel iECG-LAT maps was observed during sinus rhythm. This improved performance is attributed to the incorporation of the effects of the His-Purkinje system on the initiation of ventricular activation (Oosterhoff et al., 2016). Previously reported absolute difference for breakthrough of epicardial pacing was smaller compared to the present study (13.2–20.7 mm vs. 42.1 ± 18.6 mm) (Oosterhoff et al., 2016; Graham et al., 2019; Hohmann et al., 2019). However, previously reported absolute difference for epicardial breakthrough during sinus rhythm was higher compared to our results (75.6 ± 38.1 mm vs. 42.1 ± 18.6 mm) (Duchateau et al., 2019). Again, these differences may be explained by estimations of rhythms originating from a single ventricular focus and sinus rhythm. Thus, spatial resolution observed in this study was comparable to the earlier studies in paced complexes (Cluitmans et al., 2017; Hohmann et al., 2019). Due to the complex nature of the His-Purkinje system and the Purkinje-myocardial coupling, the implemented methods remain an approximation of the true myocardial activation and His-Purkinje physiology (Durrer et al., 1970; Myerburg, 1971; Veenstra et al., 1984).

We observed a high agreement between estimated and measured body surface potential maps, whereas the inter-map agreement was less. As the inverse problem is ill-posed, completely different ventricular activation sequences can result in similar body surface potential map waveforms, consequently we found a high agreement between body surface potential maps but a lower agreement in myocardial activation patterns.

The conduction velocities calculated on the epicardial and endocardial surfaces in this study for both the EAM-LAT maps and iECG-LAT maps were quite high (>1 m/s). However, we note that these conduction velocities are mostly determined by the velocity estimated at the surface of the myocardium. Consequently, in a Purkinje dense region, surface velocity may appear high because it also reflects the effect of the activation spread by the Purkinje fibers and not only by the myocardial tissue at the endocardial surface. Furthermore, at the epicardial surface, velocities may appear high due to the occurrence of transmural waves.

In this study, initial sites of activation were determined in the iECG method based on the observations of Durrer et al. and nine possible sites of early activation were localized (Durrer et al., 1970). Sets of these initial sites of activation were tested based on the correlation coefficient between the computed and recorded body surface potential maps, as described in more detail in the Supplementary Methods. This hypothesis was partially tested by comparing the EAM-LAT maps to the iECG-LAT maps. However, as endocardial EAM-LAT maps were often either of the right or the left endocardial surface and also did not cover the complete endocardial surface for each patient, the comparison between the number of identified EAM foci and iECG foci was hampered. This was also reflected in the absolute difference in location of identified foci of approximately 30 mm comparing iECG foci to EAM foci.

Previous versions of EDL-based methods estimating His-Purkinje mediated activation (e.g., sinus rhythm) were based on a multi-focal search algorithm over the complete endocardium and epicardium, where the first identified focus was chosen based on the highest correlation between recorded and simulated body surface potentials (van Dam et al., 2009; Oosterhoff et al., 2016). Consequently, this algorithm directly assumed that by using one focus, most of the underlying activation sequence could be ‘explained’. However, sinus rhythm, and especially narrow QRS complex sinus rhythm is an interplay between multiple activation wavefronts. Implementation of the His-Purkinje system excludes these unrealistic estimates and provides the possibility to test multiple near simultaneous foci. At the same time the initial estimation is restricted to the physiologically realistic anatomical areas and the computational burden of the iECG algorithm is minimized.

Post-processing of ECG signals, electrogram signals, and cardiac imaging influences iECG accuracy (Cluitmans et al., 2018; Graham et al., 2019). To achieve high quality EAM-LAT maps, which were used as gold standard for comparison, electrograms derived from multi-electrode catheters required re-annotation using bipolar and unipolar signals and timing to a timing reference (Cantwell et al., 2015; Graham et al., 2019). However, inhomogeneity in LAT distributions of EAM-LAT maps were observed even after re-annotation, which may have influenced the observed agreement between iECG and EAM-LAT maps. Both the epicardial and endocardial surfaces had an adequate spatial distribution of electrograms as reflected in the number of LAT per mm2 (see Figure 3C). Furthermore, the percentage of mapped surfaces was variable and some EAM procedures resulted in incomplete endocardial EAM anatomical point clouds, which affects calculated inter-map correlation coefficient (see Figures 3C,D).

Despite a considerable improvement of the iECG approach for sinus rhythm, the technique requires further adaptations and refinements that will facilitate implementation in clinical practice. Further integration of cardiovascular imaging techniques may improve performance and spatial resolution (Tung et al., 2020). Currently, the patient specific anatomical models were limited in spatial resolution by the computational models of the iECG procedure, allowing at maximum 3000 cardiac nodes, which directly affects the resolution of the cardiac anatomical model resulting in an inter-node spatial resolution of 8 ± 1mm. Diffuse or local myocardial fibrosis affects ventricular activation patterns in structurally diseased hearts. Integration of these structural abnormalities in the iECG method and refinement of the cardiac anatomical models is likely to improve imaging of electro-anatomical substrates (van Dam et al., 2009; Oosterhoff et al., 2016; Tung et al., 2020). Since electro-anatomical substrates are not limited to solely the epicardium or endocardium, iECG may allow functional imaging of such 3-dimentional substrates in patients with arrhythmias or cardiomyopathy (Tung et al., 2020). Besides diagnostic implications, non-invasive sinus rhythm iECG may play a role in the monitoring of disease progression and in sudden cardiac death risk stratification in patients with complex electroanatomical substrates, such as inherited cardiomyopathies. Eventually, reducing the number of electrodes of the body surface potential map that currently ranges from 67 to 256 electrodes, may improve clinical applicability (Hoekema et al., 1999). For EDL-based studies, also this study, the 64-electrode setup is often used (van Dam et al., 2009; Oosterhoff et al., 2016). Mathematically this setup suffices, as the number of independent signals is adequately captured using this number of electrodes and additionally, the electrodes are distributed with a high resolution in the high-gradient potential regions on the surface of the thorax (Hoekema et al., 1999).

This single center study with a small sample size included patients with structural heart disease, which may influence the generalizability of the results. Additionally, we used a set conduction velocity over the model to determine the initial estimation. This assumption may not hold in the presence of pathologies or myocardial scarring after prior ablation, but the EDL holds for homogeneous anisotropic tissue (Geselowitz, 1992).

Electro-anatomical mapping procedures and body surface potential maps were not simultaneously recorded, but in similar conditions especially concerning anti-arrhythmic drugs. During EAM, complexes were selected using dedicated Carto/Ensite EAM systems. Furthermore, sinus rhythm complexes directly following a premature ventricular complex were excluded for analysis in both the EAM and iECG-LAT map. However, a possible influence of variations in heart rate, autonomic tonus or general anesthesia cannot be excluded. The quality of gold standard EAM may have been influenced by vendor specific algorithms within the EAM systems and regional mapping by the operator during the procedure. Inherent to invasive electrophysiological studies, EAM maps consisted of electrograms recorded from consecutive sinus rhythm complexes, whereas iECG maps were derived from five sinus rhythm complexes selected from the body surface potential map.

Quantitative comparison of EDL-based iECG during sinus rhythm in patients undergoing invasive endocardial and epicardial electro-anatomical mapping showed improved agreement when compared to prior validation studies, although with considerable absolute difference in both timing and breakthrough of ventricular activation. Non-invasive iECG of both the epicardium and endocardium may prove valuable as a diagnostic tool for functional imaging of electro-anatomical substrates in sinus rhythm where activation always starts at the endocardial surface, to improve the value of the ECG in screening for cardiomyopathy and sudden cardiac death risk stratification. Future research should focus on improving accuracy and spatial resolution before implementation into clinical practice to enable imaging of functional electro-anatomical substrates.