Probabilistic Interpolation of Uncertain Local Activation Times on Human Atrial Manifolds

Objective: Local activation time (LAT) mapping of the atria is important for targeted treatment of atrial arrhythmias, but current methods do not interpolate on the atrial manifold and neglect uncertainties associated with LAT observations. In this paper, we describe novel methods to, first, quantify uncertainties in LAT arising from bipolar electrogram analysis and assignment of electrode recordings to the anatomical mesh, second, interpolate uncertain LAT measurements directly on left atrial manifolds to obtain complete probabilistic activation maps, and finally, interpolate LAT jointly across both the manifold and different S1–S2 pacing protocols. Methods: A modified center of mass approach was used to process bipolar electrograms, yielding a LAT estimate and error distribution from the electrogram morphology. An error distribution for assigning measurements to the anatomical mesh was estimated. Probabilistic LAT maps were produced by interpolating on a left atrial manifold using Gaussian Markov random fields, taking into account observation errors and characterizing LAT predictions by their mean and standard deviation. This approach was extended to interpolate across S1–S2 pacing protocols. Results: We evaluated our approach using recordings from three patients undergoing atrial ablation. Cross-validation showed consistent and accurate prediction of LAT observations both at different locations on the left atrium and for different S1–S2 intervals. Significance: Interpolation of scalar and vector fields across anatomical structures from point measurements is a challenging problem in biomedical engineering, compounded by uncertainties in measurements and meshes. New methods and approaches are required, and in this paper, we have demonstrated an effective method for probabilistic interpolation of uncertain LAT.


I. INTRODUCTION
L OCAL activation time (LAT) measurements, which indicate the onset of cardiac tissue depolarization, are part of routine clinical care for the diagnosis and treatment of atrial arrhythmias using radio frequency catheter ablation. LAT is typically determined from unipolar or bipolar electrograms recorded at different locations using an electro-anatomic mapping system, and this information is used to build a map of LAT over the atrial surface [1], [2]. This  Electrograms are recorded using a multipolar catheter which records electrical activity on the endocardial surface at one or more locations, sometimes for multiple pacing protocols and catheter positions. LAT maps can be used to calculate conduction velocity maps [3], and different pacing protocols allow the dynamic response of the electro-anatomical substrate to be inferred [4]. Although advances in catheter design have increased the number of recordings that can be made simultaneously, longer procedure times pose a risk to the patient and a cost for the health care system, which places limits on the spatial density of the mapping observations, as well as how many pacing protocols and pacing locations can be used. Data recording is therefore often sparse.
To create a LAT map, observations must be interpolated. Current interpolation techniques are predominantly linear, although a Radial Basis Function method [5] and a spline method [6] have been proposed. However, current methods do not interpolate directly on the atrial manifold, and to our knowledge no methods have been proposed to interpolate between pacing protocols. Furthermore, current methods do not properly account for observation uncertainty, nor for interpolation uncertainty. Estimating the uncertainties in conduction velocities, and therefore uncertainties in LAT maps, has been identified as an important challenge [7].
There are many potential sources of uncertainty that could influence the best estimate of LAT at a particular location. These include uncertainties in the assignment of LAT to electrograms, especially those with a complex morphology, and mismatch between catheter electrode locations and anatomical mesh for a particular recording. There have been few attempts to assign confidences to LAT observations, with some exceptions focusing on electrograms [8], [9] and uncertainty based on catheter to surface distance [10]. Besides observation uncertainty there is additional uncertainty arising from the interpolation procedure, which should be quantified so that all LAT predictions are assigned a confidence.
In this paper, we demonstrate a novel technique that accounts for these different sources of uncertainty. Our approach involves probabilistic global interpolation of heterogeneously uncertain LATs over endocardial surfaces using Gaussian Markov Random Fields (GMRFs) in a Bayesian hierarchical framework. We assign uncertainties to LAT arising from bipolar electrogram morphology by extending an existing LAT assignment method. We then estimate the uncertainties resulting from projection of measurement locations onto the left atrial surface. Using the GMRF method we construct a probabilistic LAT map over the entire left atrial surface, which accounts for heterogeneous uncertainties in LAT observations and quantifies uncertainty on LAT predictions. We extend this method to also interpolate between different pacing protocols. We apply our methods to recordings from 3 clinical cases, and demonstrate good performance using cross-validation.

A. Electrogram recordings
The recordings analyzed in this paper were part of a study involving patients with paroxysmal atrial fibrillation and undergoing first-time atrial fibrillation ablation, which is detailed elsewhere [4]. We analyzed recordings from three patients. Ethical approval was granted by the UK National Research Ethics Service (10/H0802/77) and all patients gave written informed consent for inclusion in the study. The research conformed to the principles described in the Declaration of Helsinki. None of the patients had ischemic heart disease, cardiac surgery or structural heart disease. Following femoral access and trans-septal puncture, two 8.5 French SR0 long sheaths and a PentaRay mapping catheter (Biosense Webster, CA, 1mm electrode size, 4-4-4mm spacing) were advanced into the left atrium. Decapole (St Jude Medical, MN) and pentapole (Bard Electrophysiology, MA) catheters were positioned in the coronary sinus (CS) and high right atrium (HRA), respectively. Left atrial electrophysiology was investigated using an S1-S2 programmed pacing protocol [11] consisting of a 2-beat drivetrain with a cycle length S1=470 ms, followed by a single premature extra stimulus S2<S1. The pacing protocol was delivered using a custom-built, institutionally-approved stimulator that reduced the S1-S2 interval continuously and without operator interference from 343 ms to 200 ms (or loss of capture) by reducing the current S1-S2 interval by 2% of its value (rounded to nearest ms). All pacing stimuli were delivered at a voltage of at least twice threshold and with a pulse width of 2 ms, from either the coronary sinus or high right atrium.
Bipolar electrograms (EGMs) were recorded at a sampling rate of 4 kHz using the 10 electrodes of the PentaRay catheter for each S1-S2 interval. The catheter was then moved to another location on the left atrial endocardium, and the process repeated for up to 15 locations (i.e. maximally 150 traces per S1-S2 interval). For each S1-S2 interval, we discarded any EGM traces not containing a discernible activation complex. The maximum number of discernible traces across all patient cases for any S1-S2 interval was 99.

B. Electrogram processing
A theoretical link exists between the maximal downslope of unipolar electrograms and depolarization of cardiac cells [12]. However, there is no theoretical link nor accepted best method for assigning LAT for bipolar electrograms, which are commonly used in clinical practice because they measure only local activity. Several LAT assignment methods for bipolar electrograms (BEGMs) have been reported that are sufficiently accurate for using LAT as fiducial markers for conduction velocity calculations [7]. However, the differences between these methods may be relatively unimportant compared to the intrinsic uncertainty in determining LAT for bipolar signals. We model the relation between the true local activation time (LAT EGM ) and the LAT assignment (t obs ) made from an observed EGM (derived using a specific method), as: where ǫ EGM is an error assumed to have a normal distribution with standard deviation σ EGM for tractability. We can be more confident of our LAT assignment for signals with clean single deflections than for noisy multiple deflections, so the magnitude of σ EGM should depend on the signal morphology in a robust and consistent way. Furthermore, the assessment of σ EGM should be meaningfully related to the measurement yielding t obs . For these reasons, in this paper we determine t obs from bipolar electrograms using the 'Centerof-Mass' (CoM) method [13], [14], because this has been demonstrated to accurately assign LATs to bipolar signals based on comparison with expert assignment of corresponding unipolar electrograms. In the original publications, the beginning and end of the signal were carefully defined, the signal was filtered, differentiated and rectified, and LAT was defined as t 50 : the time of 50% of the cumulative area under the rectified differentiated signal between the signal beginning and end. We extend this method in order to also give a consistent definition of σ EGM .
We summarize our modified Center-of-Mass method as four stages below, illustrated in Figure 1, and refer the reader to Supplementary Material for a schematic and more precise details for step (b). Where the original method used a 4-state machine algorithm to bracket the EGM activation complex in a very particular way, we use an amplitude threshold instead. This simplification is justified since t 50 is only a noisy observation; increasing the precision of this measurement does not reduce the intrinsic uncertainty about how t 50 differs from LAT EGM . Also, where the original method applied frequency filters to smooth signals, we use either a sliding-window average (step (b)) or a Gaussian Process (step (c)). (a) Regularize signal (Figure 1(a)); subtract the mean and normalize the amplitude of the signal V (t) to remove trends. No modifications to smoothness or morphology are made: where V (t) is the average of V (t) calculated for a time window of up to 350 ms following the pacing stimulus. (b) Bracket activation complex (Figure 1(b)). We smooth V raw (t), then differentiate (central difference), rectify, and normalize it, and then smooth again to obtain "smooth rect", which we used to bracket the complex by observing where the amplitude falls to 2.5% of the maximum. We discard V raw (t) outside the brackets. If bracketing failed then we did not proceed to (c) and no LAT observation was recorded.
(c) Smooth signal (Figure 1(c)). Before taking the rectified difference for area calculation in part (d), the signal must be smoothed. In one of three patients presented here, nearly half of the EGMs were partly 'clipped' due to thresholding by the recording equipment. Therefore we smoothed and    reconstructed V raw (t) to produce V smooth (t) using a Gaussian Process (GP) instead of a simpler smoothing method. We "censored" the clipped data and fitted a GP with Matérn 3/2 kernel (this choice prevents oscillatory artefacts in censored region) with a fixed nugget 0.001 (for smoothing) using the software GPy [15]. Since the clipping was not too severe, this enabled all of the clipped EGMs to be used, and improved the area calculation in (d). The effectiveness of the reconstruction was tested by artificially clipping complete signals; the reconstructed signals were a good match to the original. (d) Assign uncertain LAT (Figure 1(d)). We differentiated, rectified and normalized the smoothed (and reconstructed) signal from (c), and calculated the cumulative area under this signal. For the observation of LAT, we define: This particular choice (t 25 and 75 ) allows plausible features for LAT EGM to be within the ≈ 95% confidence interval implied by the error distribution in Eq. (1). Mostly importantly, this extension of the Centre-of-Mass method is motivated by our aim to have self-consistent uncertainty quantification in EGMs, such that the uncertainty is related to the EGM feature used for the observation of Local Activation Time. We return to this in the discussion (section IV).

C. Mesh processing
We obtained left atrial anatomies during the clinical procedures using an electro-anatomical mapping system (St Jude Velocity). To obtain a smoothed domain suitable for the GMRF method described in section II-D, we proceeded as follows. First, we smoothed each surface by applying a Poisson filter implemented in the MESHLAB library, [16]. Next, we determined the center location of the ostium of each pulmonary vein (PV) as the point presenting an abrupt increase of the PV diameters; we then truncated each PV at a distance of 3mm from the centre of the PV ostium. Finally, we generated a uniform triangulation with a characteristic edge length h = 2.3 mm. This procedure used the tools provided in the VMTK library, [17]. Figure 2 shows a smoothed anatomy obtained by applying this procedure. The coordinates at which electrograms were recorded (barycentres of the bipolar electrodes on the catheter) did not lie perfectly on the original (non-smoothed) left atrial mesh. Thus there was additional uncertainty in the location of LAT observations on the mesh, which we modeled as: In order to register the electrode positions with the mesh, we projected each measurement coordinate x onto the nearest coordinate on the smoothed mesh x * using a user-contributed MATLAB routine called point2trimesh [18]. Figure 3 shows this projection.
In order to meaningfully relate the uncertainty quantification to the measurement, we assume that larger projection distance d = |x−x * | ought to correspond to larger σ pos [10]. However, since the activation wavefront travels over the manifold, we need to relate uncertainty in time to a distance on the manifold. Therefore, the signal measured by the electrode was assumed to lie within an arc of length d × θ centered on x * (see inset of Figure 3), such that there is an associated LAT difference ∆t between the mesh and the catheter electrode given by: where we assume a conduction velocity of 350mm/s and fairly broad arc angle θ = 2π/10 in order not to underestimate σ pos . Assuming 95% confidence that LAT EGM is within ±∆t/2 of LAT mesh , we can write: For the GMRF method discussed below, these measurements were then assigned to the nearest mesh vertex x * * , but since this distance was small and usually already included within the arc length defined above we did not include additional uncertainty from this assignment. Upon user inspection, the nearest-point projection resulted in a few incorrect assignments on the mesh, mostly for observations with large projection distances. This happened when measurements were from parts of the original mesh that had been removed, or when a more accurate projection method was needed to distinguish between two different mesh locations of similar projection distance. Some immediately neighboring electrode locations had completely inconsistent (but clear) electrograms, demonstrating a problem with identifying the catheter position in the procedure. We labeled these inaccurate coordinates and ignored any corresponding LAT observations in Section II-D (see Table I).

D. Gaussian Markov Random Field method
A popular technique for spatial interpolation is Kriging (e.g. [19]), or (almost) equivalently Gaussian process regression [20]. To use this method we need to specify a function to give the covariance between LATs at any two vertices on the mesh. A 'standard' Gaussian process regression approach would define the covariance function to be a decreasing function of Euclidean distance between input (vertex) locations. This would not be appropriate here because Euclidean distance may not reflect the actual distance traveled by an electrical wave over the curved left atrial manifold.
To account for the mesh surface in our interpolation, we use Gaussian Markov random fields, proposed in [21], which are equivalent to Gaussian processes under some circumstances, but which allow for more complex domains and more efficient computation. We give an outline of the method here, and refer the reader to [21] for the full details. We implement this method in R [22], using the package R-INLA [23], [24]. This approach works by representing the Gaussian field on the atrial domain, denoted g(·), by a decomposition defined on a triangular domain/mesh: where each basis function ψ k is defined to be piecewise linear for any location x in the domain, taking the value 1 at the k-th vertex in the mesh, and 0 at any other vertex. The weight w k gives the value of g at the k-th vertex. Uncertainty about the field g is then described by uncertainty in the weights w 1 , . . . , w n , which are assumed to have a multivariate normal distribution.
We have noisy observations of a subset of the weights, i.e., LAT observations at a subset of mesh vertices. Conditional on various parameters (to be described shortly), it is then straightforward to estimate g at the remaining vertices conditional upon the noisy observations using standard properties of multivariate normal distributions. We now have to specify a mean vector and variance matrix for the weights. We use a constant prior expectation β for all weights. For the variance, in [21] it is shown how to construct a (sparse) precision matrix Q (the inverse of the variance matrix) for the weights, which is parameterised in terms of an unknown length scale parameter l, a smoothness parameter α (which we fix a priori), and a standard deviation parameter ω. The construction ensures a valid (positive-definite) precision matrix, with covariances between weights a function of distances along the manifold. Note that under this construction, unlike a 'standard' Gaussian process model, we do not have a closed form expression for the covariance between LAT at two different locations. By representing the unknown LAT field using (7), with the precision matrix Q, we can interpret this procedure as using a finite element method to obtain an approximate solution to a stochastic partial differential equation, whose exact solution is a Gaussian process with a Matérn covariance function with length scale l, smoothness parameter ν = α − d/2 (where d is dimension i.e. 2), and variance parameter ω on the atrial domain [21].
Our model linking the observations of LAT to the LAT field on the manifold (LAT mesh ) is defined by equations (1) and (4). Combining the errors from the electrograms and projection (assuming independence) gives: If we define y 1 , . . . , y N to be the t 50 values with associated standard deviations σ 1 , . . . , σ N , then in summary, we use the following hierarchical model: The linear predictor µ is the true LAT on the atrial manifold (i.e. LAT ≡ LAT mesh ) that we want to calculate a distribution for. The field g(·) has a Gaussian field distribution, with precision matrix Q constructed given the mesh vertex locations and parameters l, α, ω. We fix α = 2, corresponding to Markovian Matérn fields with smoothness ν = 1 and length-scale parameter l. The parameter β is introduced as an intercept, an average for the entire LAT field, so that g(·) can be assumed to be zero mean). The parameter τ is a global precision which acts to scale the observation precision σ −2 i . Although τ could be fixed to 1, we estimate its value during model fitting, so that our uncertainty assignments can be scaled to achieve a better model fit and for generality.
We use two modeling regimes: (i) LAT spatial interpolation for each S1-S2 interval independently, described in equations (9)-(11) and denoted SpatialOnly; and (ii) LAT interpolation over space and S1-S2 interval, denoted SpatialPacing. Regime (ii) involves interpolation across different pacing protocols so that the underlying random field becomes a function of both position and S1-S2 interval. We extend the model used in regime (i) to a separable space-time model, the autoregressive process of order 1, AR(1), as described in [25], where S1-S2 takes the place of time in order that the spatial fields at different S1-S2 values are correlated. We can view this as replacing equation (11) by where Q T is the precision matrix of an AR(1) process with correlation parameter ρ; see [23], [25] for more detail. It is important to note that this method is not simply interpolating between observation data from different S1-S2 values, but is correlating the underlying LAT field obtained for different S1-S2 intervals. This is, as far as we are aware, the first instance of interpolating across different pacing protocols. Fitting the model refers to the training of the parameters, and predicting refers to obtaining statistics for the linear predictor µ, which is synonymous with LAT ≡ LAT mesh . For numerical reasons we scaled LAT units into deci-seconds, and unscaled the results afterwards. Distances were measured in millimeters. The uncertain parameters are l, τ , β, ω, and for regime (ii) ρ. Penalized complexity priors [26] were chosen for l, τ , ω and ρ: a 1% probability that l was less than 100; a 1% probability that √ τ is above 2; a 1% probability that ω was above 10; a 80% probability of ρ being above 0.9. A vague uniform prior was chosen for β. To marginalize over the uncertain parameters l, τ , β, ω, we use the integrated nested Laplace approximation [27] implemented in R-INLA.
Our procedure is outlined in Figure 4.

E. Cross Validation
We predicted a LAT map for each S1-S2 interval considered. For the SpatialOnly model the data y i were observations from only the considered S1-S2 interval. For the SpatialPacing model the data y i were observations from 'triplets' of S1-S2 intervals; the central interval being the considered one.
For model validation, we used k-fold stratified crossvalidation implemented in the R package caret [28], such that validation folds were optimized for average y i . For the SpatialOnly model, the data (for each pacing) is split into k folds (subsets), and the model is trained on all data except one fold and used to predict the left-out data. This is repeated for  Fig. 4. Schematic of our LAT interpolation methodology. Various methods for assigning LAT to electrograms (LAT data) and assigning measurements to the (smoothed) anatomical mesh (spatial data), with associated uncertainties, could be used in place of our methods. Data from different S1-S2 intervals can be treated separately (SpatialOnly) or jointly (SpatialPacing). all folds. We investigated the SpatialOnly model for k = 10, k = 4, and k = 2 in order to determine how prediction was affected by the available number of mapping points.
To validate the SpatialPacing model, we apply the exact same k = 10 folds used for the SpatialOnly model to the central interval in the triplet of S1-S2 intervals, such that the validation observations are the same between both modeling regimes e.g. for S1-S2 intervals 329 ms, 336 ms, 343 ms, we perform cross validation for interval 336 ms using the same k = 10 folds used for the SpatialOnly model for interval 336 ms. This allows a direct comparison of the two models. Also with this cross-validation regime, we (i) vary the spacing between the triplet of S1-S2 intervals; (ii) combine data from triplets of S1-S2 intervals into the SpatialOnly model. We present our cross-validation results as prediction versus observation plots, and summarize them using the normalized Root Mean Squared Error (RMSE) for all predictions versus observations, calculated as: where range = max(y) − min(y) for data y i only for the considered S1-S2 interval. A perfect match of predictions E[µ] and noisy observations y i is not expected (or wanted), and RMSE scores can be highly influenced by data with large uncertainties for which the t 50 is far from LAT mesh . Therefore we also calculate the Independent Standard Error (ISE): Good predictions will have ≈ 95% of ISEs lie within the interval ±2. Both the RMSE and ISE scores, for each S1-S2 interval value, are shown in the legends in Figures 7 -8.

F. Comparison with alternative interpolation method
Our approach aims to propagate uncertainties in LAT observation to LAT predictions. A direct comparison with other interpolation methods is difficult for two reasons; there is no ground truth LAT map, and currently used interpolation methods do not incorporate uncertainties. We therefore compared the SpatialOnly model against a 3D Gaussian process with heteroscedastic observation variance because this method allowed us to incorporate uncertain LAT observations. Details are provided in Supplementary Material. Although the 3D Gaussian process application itself is somewhat novel, it is comparable with other LAT interpolation in that interpolation is not performed directly on the manifold itself. GP interpolation is very similar to Radial Basis Function interpolation [5] (we used a GP with RBF kernel) and Spline interpolation [6] (splines are special cases of Gaussian processes).

III. RESULTS
In what follows, we summarize the estimated LAT map by the posterior mean and standard deviation of the linear predictor µ ≡ LAT at each mesh vertex. The GMRF model training and prediction was performed on a HP Elitedesk desktop computer with Intel Core i5-7500 CPU (3.40GHz × 4) and 8 GB RAM (running Ubuntu 18.04.1 LTS). For the SpatialOnly model and SpatialPacing model respectively, optimizing with default parameter settings took around 60 s and 600 s, but with well-chosen starting values this reduced to 10 s and 60 s (e.g. parameter values from fitting to all data were used to initialize the optimization during cross validation). To further speed up cross validation we used the empirical Bayes approximation in R-INLA, and made predictions only for mesh points corresponding to the validation data.
We considered S1-S2 intervals between 343 ms and 280 ms, because data coverage became increasingly poor below 280 ms (possibly related to refractory tissue). This also allowed the same parameter priors to be used for all modeling. We present results for three patient cases, for pacing from the coronary sinus (CS). Data summaries are shown in Table I. Patients A and C both had fairly good data coverage overall, but the coverage varied with S1-S2 interval much more for Patient A. Patient B had very low data coverage. Figures 5 and 6 show E[µ] and V [µ] respectively for patient C for S1-S2 interval 336 ms. The LAT map E[µ] is typical of our results: the field is smooth, and high precision observations (smaller σ i ) influence and pin the LAT field more than low precision observations (larger σ i ). The heterogeneous observation uncertainty also influences the σ ≡ V [µ] map, demonstrating lower prediction certainty in regions of low precision observations. In the lower part of Figure 5, some observations (yellow spheres) appear based on physical considerations of the propagating wave to be poorly assigned to the mesh, and should probably have been identified beforehand and ignored. Nonetheless, the model mostly overlooks these problems due to influence from other more precise observations. The GMRF method is also able to smoothly extrapolate beyond the data, although predictions too far beyond the data locations should be treated with caution in general.

Fig. 5. LAT map mean E [µ]
for patient C for SpatialOnly model with observations y i for S1-S2 interval 336ms. The spheres show LAT observations y i ≡ t 50 . Sphere sizes increase with observation variance. Higher precision observations pin the LAT map more strongly. There appear to be some observations that were poorly assigned to the mesh (lower figure, yellow spheres), which should ideally be removed followed by refitting.
A. SpatialOnly model: single S1-S2 intervals Figure 7 shows the k = 10 cross-validation results for the SpatialOnly model, with the RMSE and ISE statistics shown in the legends against the corresponding color for the S1-S2 colorbar. The RMSE statistics for k = 10, 4, 2 folds, for both the GMRF model and GP model, are summarized for comparison as averages over S1-S2 intervals in Table II. Bar plots showing the breakdown of these results against S1-S2 interval are given in the Supplementary Material. Figure 7(a) shows good prediction versus observation for all S1-S2 intervals in one patient, whereas Figure 7(c) shows that the prediction appears to get worse for low S1-S2 values in another patient. The RMSE scores in Figure 7(a) are approximately half those in Figure 7(c), and the goodness of fit is further suggested by better ISE coverage in the ±2 range for the former.
The top three rows of Table II show that when the size of the cross-validation folds increase, the RMSE scores increase but the effect is very minor, demonstrating good robustness of the spatial interpolation method to decreasing amounts of data. However, for Patient B (the least data) both figure 7(b) and Table II demonstrate that prediction is much more difficult with sparse spatial coverage (although the ISE plots in Figure 7(b) suggest that the poorer predictions are compensated for by a  higher uncertainty). Interestingly, predictions with k = 2 folds for Patients A and C are still much better than predictions for k = 10 folds for Patient B, suggests that spatial coverage overall is more important that the absolute number of data points. This is because withholding validation data where the data is already sparse may result in bad extrapolation. The bottom three rows of Table II show the cross-validation RMSE summaries for the simpler 3D Gaussian process model, where the same folds from the GMRF spatial model crossvalidation were used. In all cases, our model performed better than the 3D GP method, although the GP model performed quite well overall when enough data was available. However, the GP model performance decreased much faster as the number of folds decreased, sometimes due to the inability to fit the data set well. However, as shown in the Supplementary Material, the predicted uncertainty in the GP model is not sensible given the physical system that is being modelled. This is because the interpolation can take no account of the manifold, and low uncertainty can be achieved where (c) Patient C cross-validation. Good consistent data coverage, better performance for higher S1-S2 intervals.
the distance to data 'through the atrium' is small, even though the connecting distance 'across the atrium' is large. Conduction velocity calculations (shown in Supplementary Material) demonstrate that the GMRF method outperforms the GP method substantially; the GP conduction velocity field is riddled with problematic artefacts and unphysical prediction of the electrical wave propagation, whereas the GMRF method yields very sensible results due to the interpolation being performed on the manifold itself. B. SpatialPacing model: multiple S1-S2 intervals Figure 8 shows the k = 10 cross-validation results for the SpatialPacing (multiple S1-S2 interval) model, where for each S1-S2 interval the exact same folds for the SpatialOnly model were used. The top four rows of Table III show the average RMSE results for the SpatialPacing model, where the triplet of S1-S2 intervals used in the model were adjusted across the four rows e.g. for predicting LAT for S1-S2 interval 310 ms, the nearest S1-S2 intervals were 304, 310, and 316 ms, the next nearest 298, 310, and 322 ms, and so on. Both Figure  8 and Table III demonstrate substantially better performance than the SpatialOnly model in all cases. Figure 8(c) shows that the decreasing performance with decreasing S1-S2 interval for Patient C all but disappears. The breakdown of these results by S1-S2 interval is given in the Supplementary Material. For the SpatialPacing model, the predictions seem to be just as good for S1-S2 intervals far apart as for S1-S2 intervals very close together, and certainly we have shown that for S1-S2 intervals of around 25 ms (a range of S1-S2 intervals of 50 ms) the SpatialPacing model works well and outperforms the SpatialOnly model. This is due to the ability to borrow information about the spatial activation pattern from other pacing protocols, and combine it appropriately in the modelling. This point is emphasized in the bottom four rows of Table  III, which show the average RMSE cross-validation results for combining different pacing data into the SpatialOnly model. In this case, the RMSE scores are higher and get worse as the S1-S2 interval becomes longer (for Patient B, the result for the next next next nearest neighbours is still significantly better than for modeling using data from only one pacing, due to the low amount of data in this case).

A. Uncertainty in LAT observations
Our overall methodology is modular, so other methods for electrogram analysis and data registration to the atrial mesh could be used with our interpolation method. As with this study, uncertainty quantification should be consistent with how Local Activation Time and catheter position were observed. We assumed a normal distribution to LAT in order to combine uncertainties from electrograms and position, and for computational tractability; other choices may be better. For LAT assignment to electrograms, we extended a 'Center-of-Mass' method [13], [14] that does not depend on unreliable fiducial markers, such that uncertainty depended on electrogram morphology consistently with how LAT was assigned. Our uncertainty definition was based on experience with the recordings used in this study, and the LAT distribution embraced a range of plausible features for the 'true' LAT. A more systematic calibration would optimize the relative consistency of electrogram uncertainty from one recording to another, and identify suitable timing thresholds.
In clinical data, electrograms can be highly fractionated such that there are multiple activation complexes that need to be separated; see [29] and [30]. In these cases, there is additional uncertainty associated with bracketing (separating) the complexes in the signal. It will be important in these cases not to overlook this 'bracketing uncertainty'; LAT assignment will depend heavily on how activation complexes are separated in highly fractionated signals, and it would be better to define large uncertainties than to be overconfident.
Mismatch between the electrode location provided by the electro-anatomic mapping system and the atrial mesh is a further source of uncertainty, which we have also included in our study. Registration of electrode locations with the anatomical mesh is a difficult problem because the anatomical map is produced at a different time to the measurements. Our method takes this uncertainty into account in meaningful way by directly relating the error in the registration method and the corresponding uncertainty in LAT that could likely occur. However, more accurate and robust methods of registration, such as [31] and [32], which deal both with improving the accuracy of registration as well estimating errors, could be used within our framework.

B. GMRF method for interpolation
The GMRF method can include uncertainty in LAT measurements and yield probabilistic LAT maps. In this framework the LAT observations are not taken to be the ground truth, but rather are noisy observations from which the ground truth must be inferred. Notably, interpolation is performed directly on the atrial manifold. For this reason, in comparison with a 3-dimensional Gaussian process interpolation (which is qualitatively similar to previous methods), the GMRF method performed better and made more plausible estimates of uncertainty. GP interpolation did not take account of the manifold, and so resulted in problematic artefacts in the conduction velocity field. Mesh smoothing and adjustment may also affect results in important ways, since the truncation of pulmonary veins (and the left atrial appendage, if this is removed) changes the connectivity of the mesh; it is potentially very important that the actual geometry is properly taken into account as done here.
Accurate assignments of measurements to the atrial mesh is very important and will have large consequences for accuracy; inaccurate/incorrect assignments cannot necessarily be compensated for by low precision. It is difficult to determine a minimum number of mapping points required to use the GMRF method; good predictions can be expected near to training data, but not necessarily in regions of far extrapolation. However, predictive uncertainty increases far away from the data, and so the true underlying LAT may still be captured by the predictive distribution. Our results suggest that the method is robust to decreasing the number of observations, and that an even coverage of accurate LAT recordings is more important than the absolute number of recording locations.
Our method should be applicable to more complex activation patterns such as those which occur during atrial arrhythmias. However, in order to accurately interpolate such patterns, there will need to be sufficient spatial coverage and measurements will need to have high precision, such that a more complex pattern with lower uncertainty is a better fit to the observations than a smoother pattern with higher uncertainty. It may be necessary to investigate the effect of assigning different hyperparameter priors, and to fix the absolute value of the observation uncertainty in the modeling, in order to fit more complex activation patterns. It should be noted that there is no enforcement of monotonicity in the predictions, so it is possible that our method could produce implausible predictions in regions where data coverage is extremely poor or in regions where the data has been poorly assigned to the anatomical mesh.
The SpatialPacing model demonstrates the first case of interpolation across different pacing protocols that we are aware of. Our results demonstrate that information can be shared between different pacing protocols in order to improve LAT map predictions. This potentially allows for more efficient data collection in procedures. For example, it will be possible to use totally different spatial locations for each pacing protocol in the procedure, or to collect data for only a subset of pacing protocols for each placement of the recording catheter. This could radically improve the spatial coverage of data obtained during procedures with multiple protocols.
LAT maps are often patchy, but our results suggest that non-smooth patterns resulting from overlooking uncertainties should be treated with caution. The smoothness of the LAT maps produced by our method is somewhat a result of the observation uncertainties. From a physical point of view, the activation pattern should be assumed to be smooth, unless the data (with observation uncertainty included) strongly suggest otherwise. For this reason, our method of interpolation could be thought of as assuming a smooth prior for the activation pattern, such that more complex activation patterns must be enforced by the data. Probabilistic interpolation of uncertain Local Activation Times should therefore be an important new tool for constructing LAT maps for clinical purposes.

V. CONCLUSION
In this paper we have described and evaluated a novel method to interpolate uncertain Local Activation Times directly on the left atrial manifold and between pacing protocols. Our approach is to first quantify uncertainty in assigning Local Activation Times (LAT) to bipolar electrograms, as well as uncertainty in registration of catheter positions to patient specific anatomical meshes, and then to explicitly incorporate this uncertainty in LAT into a probabilistic interpolation scheme using Gaussian Markov Random Fields (GMRFs) in a Bayesian hierarchical framework. Our methods yield probabilistic LAT maps with predictive mean and standard deviation which account for heterogeneous observation uncertainty and interpolation uncertainty. We used cross-validation to demonstrate performance, showing that our method outperforms a 3D interpolation method which is comparable to existing methodologies, and demonstrated that our method can sensibly combine data from different pacing protocols to improve predictive performance. This approach enables an approach for building maps of LAT in individual patients that is robust in the face of measurement uncertainties.