{'fraction_non_alphanumeric': 0.02950940612320177, 'fraction_numerical': 0.008510301944459082, 'mean_word_length': 4.854542650007713, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gaussian Process Manifold Interpolation for Probabilistic Atrial Activation Maps and Uncertain Conduction Velocity

In patients with atrial fibrillation, local activation time (LAT) maps are routinely used for characterising patient pathophysiology. The gradient of LAT maps can be used to calculate conduction velocity (CV), which directly relates to material conductivity and may provide an important measure of atrial substrate properties. Including uncertainty in CV calculations would help with interpreting the reliability of these measurements. Here, we build upon a recent insight into reduced-rank Gaussian processes (GP) to perform probabilistic interpolation of uncertain LAT directly on human atrial manifolds. Our Gaussian Process Manifold Interpolation (GPMI) method accounts for the topology of the atria, and allows for calculation of statistics for predicted CV. We demonstrate our method on two clinical cases, and perform validation against a simulated ground truth. CV uncertainty depends on data density, wave propagation direction, and CV magnitude. GPMI is suitable for probabilistic interpolation of other uncertain quantities on non-Euclidean manifolds.


Background
Electrical activation of cardiac tissue acts to initiate and synchronise mechanical contraction. The spread of the activation wave is characterised by local activation times (LAT) and conduction velocity (CV), where CV describes the local speed and direction of the wavefront [1]. CV is a clinically important material property of cardiac tissue [2,3], with decreases in CV often associated with abnormalities such as tissue fibrosis. CV can also be used to calibrate patient-specific models of electrophysiology [4], which have the potential to be used for planning clinical interventions [5]. CV can be derived from measurements of LAT at a set of locations. In the clinical setting, LAT can be obtained from electrograms recorded using catheters placed inside the heart. However, reliable CV estimation based on uncertain estimates of LAT from noisy electrograms is difficult [1,6]. Wave collision and lines of block also present challenges for CV estimation [7,8].
In a previous study of electrical activation in the human left atrium, we identified uncertainties associated with identifying LAT from individual electrograms, as well as uncertainty associated with the registration of the recording electrode to an anatomical mesh [9]. We demonstrated that it is possible to create probabilistic LAT maps over the left atrium, which is modelled as a 2D manifold in 3D space represented as a triangle mesh. These LAT maps take account of uncertainties in both LAT measurements and the interpolation procedure itself. Importantly, the interpolation was performed directly on the manifold geometry, taking account of the physical distances travelled by the wavefront and the connectivity between points on the mesh. We used a stochastic partial differential equations (SPDE) approach linking Gaussian Markov random fields (GMRF) to a class of Gaussian processes (GP) with Matérn covariance functions [10]. This enabled us to construct random fields on non-Euclidean domains such as the manifold describing the left atrium.
Given a probabilistic LAT map on an atrial mesh, a CV map can easily be calculated from the LAT map posterior mean. A natural approach to compute how our uncertainty about LAT affects our estimate of CV is to use stochastic sampling. We can simulate samples from our posterior distribution for LAT, calculating CV for each of them, thus estimating a distribution for CV at every mesh location. However, doing this with the GMRF model used in [9] shows that this model is not suitable for estimating uncertainty about CV, because random samples of LAT are insufficiently smooth. In practical terms this leads to samples of CV that are not consistent with observations of LAT. In technical terms, the Matérn class of covariance functions depends upon a smoothness parameter ν, such that process samples are k-times mean-squared differentiable if and only if ν > k [11]. For the SPDE-GMRF, ν is constrained to be one, so that posterior LAT samples are not differentiable. Consequently, they are too rough to represent a spreading electrical wavefront, meaning the approach cannot be use to calculate the CV distribution. In this paper, we describe a new GP approach that does allow us to estimate CV distributions from a probabilistic interpolation of noisy LAT measurements.
The reason for using the SPDE-GMRF approach was that an 'ordinary GP', for which the covariance matrix is a function of the Euclidean distance matrix between 3D spatial locations, would not account for the geometry or topology of the atrial manifold. It is not possible to simply replace Euclidean distances with for example geodesic distances because the resulting covariance matrix is not positive semi-definite. However, the meaningful distance between locations for this physical system is geodesic distance on the manifold, not distance through the manifold that ignores the topology. A recent advance in reduced-rank Gaussian processes showed that it is possible to expand a stationary covariance kernel in terms of eigenfunctions of the Laplacian operator on the spatial domain [12]. The contribution of this paper is to address the problem of obtaining distributions of CV from probabilistic interpolation of LAT. We apply this recent advance to a 2D manifold in 3D space, so we are able to form a Gaussian process on the manifold using any stationary covariance function, including those which lead to differentiable samples. This allows us to calculate CV distributions by probabilistic interpolation of LAT. This paper is organised as follows. We first explain how to interpolate LAT using ordinary Gaussian processes. We then review how to expand the covariance function as in [12], but applied to a 2D triangulated manifold representing the atrium, such that we can simply replace the covariance matrix in ordinary GP regression to obtain a GP confined to the manifold. Using clinical data, we perform LAT interpolation and obtain corresponding probabilistic CV maps. We also validate our methodology against a simulated ground truth. We find that CV uncertainty depends on LAT observation uncertainty, data density, alignment of observations with respect to wave propagation direction, and on CV magnitude. We then discuss future work made possible by solving the geometrical aspect of probabilistic interpolation on atrial manifolds.

Data collection
As proof of concept, we present two clinical cases in this paper. Additional details on data collection are given in the Supplementary Material. For the first case, Patient 1, (identical to Patient A in [9]) observations were collected using an S1-S2 pacing interval: a 2-beat drive train with a cycle length of S1 470ms followed by a single premature extra stimulus S2, where S2 ¡ S1. For the second case, Patient 2, a constant pacing at cycle length 500ms was used, i.e., S1 only. For Patient 1, electrograms were processed to obtain uncertain LAT as described in [9]. For Patient 2, visual inspection of each point was undertaken by one of the authors to assign LAT to the earliest sharp deflection, without an assigned uncertainty based on the electrogram analysis.
In both cases, a triangulated mesh of the left atrium, obtained using the respective mapping systems, was processed to clip away the pulmonary veins and mitral valves to leave holes, and the triangulation was smoothed (using mainly surface-preserving Laplacian smoothing, element sub-division, and short-edge collapsing) to be reasonably regular, with ∼ 10 4 non-intersecting mesh faces. This yields an interpolation resolution mesh, in which LAT observations were assigned to the nearest vertex, and given an additional observation uncertainty arising from this mesh assignment as in [9]. Note that our method does not require observation-specific uncertainty for LAT, or specific methods for determining LAT and assigning observations to a mesh.

Gaussian process regression
A function f , for which only noisy observations y i are available at spatial locations x i , can be modelled by a Gaussian process as follows: where i is the (heteroscedastic) noise for each observation. The covariance function or "kernel" k specifies the covariance cov(f (x), f (x )) := τ 2 k(x, x ) as a function of location. Here, the mean of the Gaussian process is set to zero. We can generalize by using a non-zero parametric mean function, but by centering and scaling the data y to be mean zero, variance one, we find that a zero mean is sufficient and simplifies the exposition. A GP is 'fit' to data by conditioning upon the observations D = {x i , y i } n i=1 and optimizing τ 2 and any kernel hyperparameters to maximize the log-likelihood. Here σ i is obtained by pre-processing LAT and mesh assignment rather than by optimization, although such optimization would be possible.
Posterior predictions (i.e., using the GP conditioned upon the data) are multivariate Gaussian distributions specified by the posterior mean and variance. At input (e.g., location) x * , the posterior mean and variance is given by: where K is an n × n matrix containing the covariance between all pairs of observations τ 2 k(x i , x j ); y is the vector of observations y i ; and Σ ≡ diag(σ 2 1 + η, . . . , σ 2 n + η) contains the observation error variances σ 2 i (for homoscedastic noise, all σ i are equal) plus a 'nugget' η representing additional observation uncertainty common to all observations (η is optimized with other hyperparameters). The cross-covariance between observations and predictions is given by the vector k * , and τ 2 k(x * , x * ) is the covariance between prediction locations (both are matrices when x * contains multiple locations). By predicting at multiple locations simultaneously, e.g., all points on the atrial mesh, we obtain expressions for the posterior mean and variance at all locations.

LAT interpolation with a Gaussian process
We assume the underlying LAT field is a smoothly varying function of location and model it as a Gaussian process: LAT ≡ f (x) ∼ GP(0, τ 2 k(x, x )). We model the relationship between the noisy LAT observations y and the 'true' LAT as y i = f (x i ) + i where i ∼ N (0, σ 2 i ), as in Eq (1). LAT observation y i is assigned to a coordinate on the atrial mesh x i .
For the GP kernel, we choose the Matérn 3/2 kernel (ν = 3/2), which gives a once (mean squared) differentiable process. This choice was made because CV is obtained from the derivative of LAT, but we do not want to make assumptions about the smoothness of CV. The Matérn 3/2 kernel is a stationary kernel that, when working in Euclidean space, depends only on the Euclidean distance r ij = ||x i − x j || 2 between observations i and j, i.e. k(x i , x j ) ≡ k(r ij ), and can be written, with length-scale hyperparameter l, as: In [9], we jointly interpolated LAT across both space and S1-S2 pacing using an order one auto-regressive model (i.e., AR(1)) for the dependence on S2. We do not interpolate over S2 in this paper, but could do so using a separable covariance structure such as k spatial-pacing where k pacing (S2, S2 ) is a suitable standard kernel (e.g. squared exponential), simplifying the approach from [9].

Gaussian processes on a manifold
The standard GP approach described above assumes, due to the kernel, that correlations are a function of Euclidean distance. Since electrical activation advances over the atrial manifold, instead of working with Euclidean distances ||x − x || 2 in the domain R 3 , we want to limit the domain to locations on the atria, which we denote by Ω, and assume correlations are a function of geodesic path length on that domain. Unfortunately, it is not possible to simply replace the Euclidean distance in a kernel such as Eq. (3) with the geodesic distance and still form a valid covariance matrix, as this leads to non-positive definite kernels [13].
Here, we take advantage of a recently developed [12,14] reduced-rank method for expressing a stationary kernel in terms of the kernel spectral density and eigensolutions of the Laplacian operator on Ω, that allows us to bypass the need to explicitly define a covariance function. We choose a spectral density, and then define a stochastic process with that spectral density, which then implicitly defines a covariance function on Ω (by the Wiener-Khinchin theorem).
We use a finite dimensional model of the form where the f k are parameters to be learnt.f (x) is a (finite dimensional) Gaussian process with covariance function given by the basis function expansion The basis vectors φ k (x) are chosen to be eigenfunctions corresponding to eigenvalues λ k of the negative Laplacian operator on the domain Ω, i.e., {λ k , φ k (·)} are solutions of the system with 0 ≤ λ 1 < λ 2 < . . .. Note that this basis is independent of our choice of covariance structure/spectral density, allowing us to try different structures without needing to recompute the basis vectors. We use Neumann boundary conditions specifying zero derivative normal to the mesh edges (pulmonary veins and mitral valve), as these do not specify LAT in advance of the interpolation (in contrast [12,14] used Dirichlet boundary conditions). S(w) is a spectral density [11] which we are free to choose; it is S that determines the nature of the covariance structure of our model. In this paper, we use the spectral density corresponding to the Matérn family of covariance functions, which is given by where ν is smoothness, l is length-scale, D is dimensionality (D = 2 for a 2D manifold in 3D space), and Γ is the Gamma function. The spectral density of our modelf is only approximately S(w), but converges to S(w) as M and Ω grow [12]. The motivation in [12] for using expansion (4) is to reduce the computational complexity of hyperparameter estimation, as the basis is independent of the covariance function. In contrast, our motivation for using this approach is that the solution of Eq. (6) can be computed numerically on a discrete mesh, even for non-Euclidean domains, allowing us to define a GP model directly on the manifold defined by Ω. In particular, we have chosen to model LAT on an atrial manifold using a stochastic process with Matérn spectral density.

Eigenfunction calculations
To apply this kernel expansion (5) to a 2D manifold in 3D space, we represent the manifold as a discrete triangle mesh with N vertices and represent the Laplacian operator with the 'mesh/cotan Laplacian' (see [15]) given by: where s represents an arbitrary scalar field defined at vertices on the mesh. Figure 1 (left) represents the terms in eq. (8). Applied to a triangular mesh, the N × N operator matrix L represents the Laplacian with Neumann boundary conditions at the mesh edges (mitral valve and pulmonary veins). We calculate the Laplacian matrix L neglecting the areas A i , and solve the eigenvalue problem Lφ = λAφ where A is a diagonal matrix of areas A i (see [16]), solving for the M smallest eigenvalues and corresponding eigenvectors using the eigsh function from the Python package Scipy [17]. We obtain M eigenvectors such that φ k (x i ) ≡ φ i k , where k = 1 . . . M index the eigenfunction and i = 1 . . . N index the mesh vertices. These eigenfunctions belong to the spatial domain represented by the mesh, and as such the topology (connectivity) and geometry (distances) of the atrial manifold are properly accounted for. Note that the eigenproblem only needs to be solved once for a particular mesh. We can then construct all of the covariance matrices needed for GP regression. We refer to this method as Gaussian Process Manifold Interpolation (GPMI). Our code for applying GPMI to atrial manifolds can be found at [18].
To represent the atrial mesh smoothly, we use an interpolation resolution mesh of ∼ 10 4 vertices which corresponds to an edge length ∼ 1mm. In order to ensure a fine discretisation of the Laplacian, and also to allow the calculation of eigenfunction gradients (required for the posterior of the GP gradient; see Section 22.6), we subdivide the interpolation mesh: each triangular element (3 vertices) is subdivided into 9 elements by creating two additional vertices along each edge and one additional vertex at the element centroid, as shown in figure 1 (right). The Laplacian and the corresponding eigenfunctions are calculated for this finer mesh.
For each original mesh face, there are now 10 vertices which lie in the plane of this face. To calculate the gradients of the eigenfunctions at a face centroid, we first rotate these (finer mesh) vertex coordinates into a 2D frame, and then fit a cubic polynomial to the vertex eigenfunction values and 2D coordinates using least squares regression. We then evaluate the gradient of the eigenfunctions at the centroid using the derivative of the polynomial expression, and rotate this (vector) gradient back into the original 3D coordinate frame. This allows us to calculate the eigenfunction gradients at the centroids, even when the eigenfunctions vary non-linearly across a face. For simplicity, we retain the eigenfunction values at the mesh vertices and face centroids of the original mesh only (i.e. we discard the vertices from subdividing the face edges), which allows us to assign LAT observations and make LAT predictions at both vertices and centroids.
Boundary conditions at the edges of the mesh (holes representing the mitral valve and pulmonary veins) will leave unwanted artefacts in the interpolation. To deal with this problem, after forming the interpolation resolution mesh and prior to subdivision as described above, we append tubes of triangular elements to the atrial holes, consisting of 10-20 layers of elements (depending on the mesh resolution) and extending away from the mesh centre of mass, using a custom algorithm described in the Supplementary Material. This takes the boundaries away from the 'real atrium' and removes the unwanted boundary effects on the interpolation. Note that the topology is unchanged by this procedure. After subdivision as described above, we retain only the eigenvector entries corresponding to the original mesh vertices and centroids before extension. A figure of the mesh extension result for Patient 2 is given in the Supplementary Material.

Conduction Velocity Distribution
Conduction velocity is the inverse of the LAT gradient. The distribution for either the LAT gradient or CV could be calculated by first generating samples from the posterior distribution for LAT, and then calculating the desired quantities for each sample. However, since applying a linear operator (such as the gradient) to a Gaussian process results in another Gaussian process, the gradient of LAT is also a Gaussian process. Therefore, we can instead sample directly from this gradient distribution. The posterior distribution for the gradient of LAT is given by similar expressions to the posterior distribution for LAT (eq. 2), except that the covariance entries involving prediction locations require derivatives of the covariance function with respect to these prediction locations: When the kernel is expressed as in equation (5), the necessary kernel derivatives are given as: The derivatives of the eigenfunctions can all be calculated in advance of the interpolation, after calculating the eigenfunctions themselves (see Section 2.5). Note that it is therefore possible to include uncertain observations of the gradient of LAT into the modelling. Given the gradients of the eigenfunctions at the face centroids (see Section 2.5), the posterior distribution of the gradient of LAT can be calculated at all centroids. However, since neither magnitude or inverse are linear operations, we calculate the distributions for the magnitude of the gradient of LAT and CV by sampling from the posterior LAT gradient and applying the operations to the samples. At each centroid, we generate 2000 samples from the multivariate normal distribution for the components of the gradient of LAT given by equation (9). We calculate the standard deviation, 9th, 25th, 50th, 75th, 91st percentiles for the magnitude of the posterior samples of the gradient of LAT (the percentiles can be inverted to obtain CV percentiles). We use the posterior mean of (9) to represent prediction of the magnitude, and statistics from sampling for prediction uncertainty.

Simulations
To quantitatively evaluate our methodology, we compare predictions of LAT and CV against a simulation 'ground truth'. A computational mesh with a regular triangulation with edge length of 215µm was obtained for Patient 1 (1047446 vertices). Simulations were performed with the Cardiac Arrhythmia Research Package (CARP), an electrophysiology solver suitable for high-performance computing [19]. We simulate S1 pacing from the coronary sinus for S1 750ms, using the modified Mitchell-Schaeffer ionic model with heterogeneous parameters to generate a spatially heterogeneous conduction velocity field (details in Supplementary Material). To create a 'virtual patient' from these simulations, we downsample the simulation mesh to an interpolation resolution mesh, and match the interpolation mesh vertices to the nearest simulation mesh vertices. For each interpolation mesh vertex, we record LAT from the matching simulation vertex.
We calculate CV for the simulation in two ways. Firstly, we calculate CV element-wise at each interpolation centroid from LAT at the face vertices ('element' method): given a face defined by anti-clockwise vertices {i, j, k}, face area A, face normal N , edge vectors e ab , with LAT at vertex l given by t l , the centroid LAT gradient and CV are given by (× denotes the vector product) Secondly, since the CV field is heterogeneous and CV often varies on a length-scale smaller than the distance between observations, we also calculate CV taking account of a larger area ('wave' method). For the 20 vertices nearest to where CV is to be calculated, the vertex coordinates are projected into 2D such that the geodesic distances are optimally conserved. A least squares linear fit to LAT is used to calculate a gradient corresponding to plane wave propagation [8]. We do this for 1000 random centroids locations, optimized to be well spaced over the mesh using a a maximin distance criterion over 10, 000 designs. Note that 20 vertices corresponds to at least two 'rings' of vertices around a centroid.
To validate our approach for predicting LAT and CV, we calculate the normalized root mean square error (nRMSE), reported as a percentage, and independent standard error (ISE) for all predictions E[z i ] (with variance V[z i ]) versus ground truths z i , where nRMSE (%) and ISE are given by: where range = max(z i ) − min(z i ).
For LAT, about 95% of ISE scores should lie within the interval ±2. We emphasise that LAT, and therefore (implicitly) the gradient of LAT, is what is directly modelled in our framework. Therefore, for comparison with the ground truth, we predict and compare against the magnitude of the gradient of LAT, rather than against CV, using the standard deviation of the magnitudes of LAT gradient samples. This prevents the nRMSE and ISE scores from being heavily distorted by the inversion involved in calculating CV from the LAT gradient (see Section 4). We calculate these scores for the LAT gradient for comparison with the two methods of calculating CV outlined above. The meaning of the terms in equation (12) for the magnitude of the LAT gradient must be explained: we use the magnitude of the posterior mean of the LAT gradient as E[z i ], and the standard deviation of the magnitude of the posterior LAT gradient samples as

Results
The number of interpolation mesh vertices and centroids, and the number of observations, are given in Table  1. For the virtual patient, the number of observations will be given in the text for each case. We use M = 256 basis functions, which we determined to be sufficient for modelling LAT (see Section 3 in the Supplementary Material).  Validation results for the virtual patient are given in Table 2, for prediction of LAT and the magnitude of LAT gradient. For each table row, LAT observations were chosen by picking n random vertices and corresponding LAT values, and corrupting these values with additive random noise i ∼ N (0, 1). The nRMSE and ISE scores were calculated and averaged across 10 different observation designs for each n. Figure 2 shows LAT interpolation for Patient 1, overlaid with CV vectors coloured by the LAT posterior mean (top) and standard deviation (bottom) at the face centroids. The LAT observations y i and noise standard deviations σ i are shown as colored spheres. LAT prediction uncertainty depends on the distance from the observations and on the observation uncertainty. These interpolation results are comparable with those in [9]. Figure 3 shows the LAT interpolation for Patient 2. In this case, the observation noise only arises from assigning observations to mesh vertices and not from LAT assignment to electrograms, resulting in a very small observation noise. Prediction uncertainty is significantly lower than for Patient 1 due to the increased spatial data coverage.

LAT interpolation
The validation scores for LAT prediction for the virtual patient show that ISE coverage (|ISE| ≤ 2) is consistently ∼ 95%, and the nRMSE is already less than 1.0% for 250 observations. Note that these scores were comparing all LAT predictions against all ground truth LAT values. These scores demonstrate excellent prediction of LAT using the GPMI method. Figure 4 shows CV for the virtual patient. In the top row, CV is calculated at centroids using the 'element' method equation (2.7) with known (ground truth) LAT at the mesh vertices. The second and third rows show CV prediction using 1000 random observations of LAT at mesh vertices, with random noise ∼ N (0, 1) corrupting these observations for the third row. The predictions capture the heterogeneous CV field produced in the simulation, though the features are more smoothed in the prediction, especially where the true CV varied on a lengthscale shorter than the resolution given by the distances between observations. The predicted CV field is slightly smoother in the case with noisy LAT observations. Overall, the predictions seem to perform very well, capturing even the locations of wave collision and stimulus location. Figure 5 shows the prediction uncertainty (interquartile range) for CV magnitude, corresponding to the predictions in Figure 4. There is significantly more uncertainty in CV for prediction using noisy LAT observations (the colorbar scales are different). Higher CV correlates with higher prediction uncertainty (as is also the trend for Patient 2, discussed below). Table 2 gives the validation scores for the magnitude of the LAT gradient for the virtual patient. The scores are very similar for both methods of calculating the ground truth. Notably, the nRMSE scores are an order of magnitude larger for |∇LAT| prediction than for LAT prediction, although the ISE scores suggest that the truth is captured well by the uncertainty. For 250 observations, the nRMSE score for |∇LAT| is less than 10% (and the nRMSE score for LAT is less than 1%. We propose that 250 observations (spaced over the atria) could be considered a lower bound for having confidence in CV predictions calculated with our methodology. Figure 6 shows prediction vs observation plots (left) and ISE plots (right) for the virtual patient for 1000 noisy observations, where the ground truth is for |∇LAT| calculated by the 'wave'. Points are ordered so that larger standard deviations are plotted in the foreground. The error bars in the prediction vs observation plots correspond to the 9th and 91st percentiles. The uncertainty mostly captures the ground truth. Notably, the uncertainty is highest where |∇LAT| is overestimated (i.e. CV is underestimated) relative to the ground truth. This is likely because regions of fast CV cannot be captured when uncertainty is high.

Conduction velocity
For Patient 1, there are only 61 LAT observations, much less than the 250 minimal number of observations indicated from the simulation studies. These observations are also highly clustered such that large sections of the mesh have no observations at all (see Figure 2). Therefore, we do not present plots of CV magnitude for Patient 1. Figure 7 shows CV calculated from the posterior of the LAT gradient for Patient 2, coloured by the CV magnitude (top) and the interquartile range (bottom). The colored spheres represent these values (interpolated from neighbouring faces) at observation vertices. The predicted CV magnitude is sometimes very high far from the observations, where the GP flattens out beyond the data (shallow LAT gradients correspond to high CV), and in regions of wave collision where LAT is interpolated across two different wavefronts (the calculated value of CV does not represent wave propagation).
CV uncertainty depends highly on distance from observations, but also on the direction of the (predicted) CV field, so that physiological and more certain CV predictions are made where there are several nearby observations aligned along the direction of (predicted) wave propagation. This makes physical sense as observations aligned along the direction of changing LAT are required to inform the gradient of LAT, and is especially useful since we do not know the direction of propagation in advance of the interpolation.

Discussion
In this study we used Gaussian processes to interpolate LAT directly on atrial manifolds in order to calculate CV distributions. The results show that point-wise uncertainty in LAT interpolation is not sufficient for inferring uncertainty in CV. Rather, CV uncertainty depends on LAT observation uncertainty, data density, and the presence of sufficient observations along the direction of wave propagation. Also, larger CV will generally have higher associated uncertainties due to the inverse relationship between the gradient of LAT and CV, such that the same uncertainty in LAT will have much larger effects on CV for shallow LAT gradients (for a given spacing of percentiles for the magnitude of the gradient of LAT, the percentiles on the distribution for CV will become further apart as the gradient of LAT is decreased). Precisely measuring fast CV is therefore difficult. Uncertainty in CV restitution curves could be important for calibrating atrial models, e.g., in [20,4] where the CV uncertainty could weight the comparison between simulated and observed restitution curves.
We have focused here on solving the geometry aspect of the problem of probabilistic interpolation on atrial manifolds, which had only been possible previously for the special case of Gaussian Markov random fields (as in our previous publication [9]). The method presented here, Gaussian Process Manifold Interpolation (GPMI), allows us to solve the problem of calculating CV statistics from a probabilistic LAT interpolation. It would also be suitable for probabilistic interpolation of other quantities such as voltage and tissue parameters. The flexible modeling framework presented here means that observations of LAT gradient could also be included into the modeling. This could allow both LAT and LAT gradient observations (for example, from multi-polar catheters where consideration of wave propagation is used to obtain precise CV values [8,6]), to be used in a unified and consistent modelling framework. It may also be possible to include linear constraints on interpolation via the covariance function [21], and to directly include physics using latent force models as in [22].
In theory, it should be possible to apply this modelling framework to a 3D domain such as the ventricles. Besides the increase in the number of basis functions that may be required to accurately expand the kernel [12], the main difficult may be in the boundary conditions at the surface of the domain. Some form of mesh extension analogous to that done here may be possible, or perhaps different boundary conditions might be applied. One possibility may be to utilize a separable kernel structure, similar to the 'spatial-pacing' model mentioned in Section 22.3, whereby the 'thickness' direction is accounted for with a separate kernel that multiplies the manifold kernel.

Limitations
LAT interpolation is difficult in part because of genuine discontinuities where wavefronts collide and where lines of block may exist, and also because of the quasi-monotonic nature of the underlying activation field. These difficulties are not solved with a Gaussian process, which may further introduce artefacts such as reverting to the mean in regions of extrapolation and oscillations in larger gaps between observations. Measurements of CV made by fitting plane waves to LAT observations can explicitly model a propagating wavefront, but the LAT interpolation here does not. Furthermore, CV calculations using LATs measured from individual placements of a multipolar catheter have the advantage that the distance between recording sites is fixed and known, whereas LAT interpolation and CV prediction from combining observations for many separate placements of a catheter (as done here) is potentially subject to large inaccuracies, i.e., some non-physiological CV predictions simply reflect the data, and it may be hard to correct for this without making strong prior assumptions about CV.
Accurately predicting the gradient of LAT is more challenging than predicting LAT, as noisy observations of LAT provide limited gradient information, particularly when sparsely spaced. Furthermore, the magnitude of the LAT gradient is not a Gaussian process. Nonetheless, we find that the calibration of the predictions (the coverage achieved by credible intervals averaged over all predictions) is reasonable for the posterior distribution of the LAT gradient magnitude.

Conclusion
We have introduced the Gaussian Process Manifold Interpolation (GPMI) method which allows us to perform probabilistic interpolation directly on atrial manifolds using a reduced-rank Gaussian process that expresses the covariance function in terms of eigenfunctions of the Laplacian operator on a mesh triangulation of the manifold. This allows us to interpolate uncertain local activation times with a differentiable model, which then allows us to calculate probabilistic estimates of conduction velocity maps. In general, this method allows Gaussian processes to be used to model data on irregular spatial domains such as the left atrium, while directly accounting for the topology of the domain and the physically meaningful distances that correlations ought to depend upon.
Ethics: Clinical data were recorded as part of routine care for first-time ablation for atrial fibrillation. Ethical approval was granted by the National Research Ethics Service (10/H0802/77) for Patient 1, and (REC number 18/HRA/0083) for Patient 2.
Contributions: SC, RW, JO and RC conceived the study. SC designed and implemented the main methods. CC, CHR and DOH processed the raw clinical data. SW and MN performed the clinical data collection. All authors read and approved the manuscript.   Figure 4. Note that the color scales are different for each case, due to differences in prediction uncertainty for CV magnitude.

Mesh extension
In order to remove the effects of the Laplacian boundary conditions from the interpolation on the atrial mesh, we extend the mesh with elements as in Figure 1. This does not change the mesh topology, and the shortest distance between points on the original mesh is almost completely unchanged (with small exceptions caused by non-planer clipping of the clinical mesh to produce the holes in the interpolation mesh; for planar or convex edges, the shortest distances would be completely unchanged). The algorithm is as follows. For each hole in the mesh (e.g. mitral valve), the vertices belonging to the 'edge' describing this hole are identified (face edges at the mesh edge belong to one face only; this can be used to identify vertices on the mesh edges). These vertex coordinates are used to fit a plane, and the plane normal n sign is adjusted such that n points away from the mesh centre of mass. Vector n points in the direction that layers of elements will be added. Then, starting on any vertex in the edge and proceeding either clockwise or anti-clockwise until all edge vertices have been considered, the following is done: 1. calculate the centre x p between current vertex x i and neighbouring vertex x j 2. add a new vertex at position x n = x p + n × |x j − x i | tan(60 deg)/2 3. append this new vertex x n to the list of all vertices 4. append a new element connecting x i , x j , x n to the list of all elements Now the edge will look like the teeth of a saw. New elements must be added to the list of all elements that connect the newly created vertices to each other, and to the original vertex nearest both of them. This concludes the creation of a single layer of new elements. This process can be repeated to build up sufficient layers to prevent the boundary conditions from producing artefacts in the interpolation (the plane normal vectors can be recalculated for each layer, or the original vector used for all layers).

Explained Variance
The percentage of variance captured by the first m basis functions can be calculated as where S is the kernel spectral density evaluated at the square root of the eigenvalue λ j (kernel hyperparameters have been optimized). Figure 2 shows, for the model fit of the virtual patient with 1000 noisy observations and M = 350 basis functions, that the first 256 eigenfunctions capture 98.8% of the variance. Therefore 256 eigenfunctions are more than sufficient for LAT interpolation.