Space–time shape uncertainties in the forward and inverse problem of electrocardiography

Abstract In electrocardiography, the “classic” inverse problem is the reconstruction of electric potentials at a surface enclosing the heart from remote recordings at the body surface and an accurate description of the anatomy. The latter being affected by noise and obtained with limited resolution due to clinical constraints, a possibly large uncertainty may be perpetuated in the inverse reconstruction. The purpose of this work is to study the effect of shape uncertainty on the forward and the inverse problem of electrocardiography. To this aim, the problem is first recast into a boundary integral formulation and then discretised with a collocation method to achieve high convergence rates and a fast time to solution. The shape uncertainty of the domain is represented by a random deformation field defined on a reference configuration. We propose a periodic‐in‐time covariance kernel for the random field and approximate the Karhunen–Loève expansion using low‐rank techniques for fast sampling. The space–time uncertainty in the expected potential and its variance is evaluated with an anisotropic sparse quadrature approach and validated by a quasi‐Monte Carlo method. We present several numerical experiments on a simplified but physiologically grounded two‐dimensional geometry to illustrate the validity of the approach. The tested parametric dimension ranged from 100 up to 600. For the forward problem, the sparse quadrature is very effective. In the inverse problem, the sparse quadrature and the quasi‐Monte Carlo method perform as expected, except for the total variation regularisation, where convergence is limited by lack of regularity. We finally investigate an H1/2 regularisation, which naturally stems from the boundary integral formulation, and compare it to more classical approaches.


| INTRODUCTION
Electrocardiographic recordings at the body surface, such as the standard 12-lead electrocardiogram (ECG), are a direct consequence of the electric activity of the heart. The spatial resolution of such recordings depends on the number of electrodes placed on the chest, ranging from no more than 10 of the standard ECG to hundreds of electrodes composing high-density body surface potential maps (BSPMs). Along with an accurate description of the torso anatomy, BSPMs are sufficiently informative to enable a non-invasive characterisation of cardiac electrophysiology, termed ECG imaging. 1 ECG imaging technologies have been extensively validated in experimental, animal and, more recently, clinical settings, with promising results. 2 The ECG imaging problem can mathematically be formulated as an inverse problem. The most classical formulation of it is associated with a potential-based forward problem, which amounts to determining the BSPM on the chest from the electric potential at a surface enclosing the heart, for example, the epicardium. The inverse problem of electrocardiography consists therefore in recovering the epicardial potential from the BSPM. 3 As typical for inverse problems, however, the ECG imaging problem is severely ill-posed in the sense of Hadamard, since arbitrarily small perturbations of the BSPM, such as noise, may yield large variations in the epicardial potential. It thus requires regularisation to penalise non-physical solutions and a strategy to optimally select the associated regularisation parameter. For the ECG inverse problem, several kinds of regularisations have been proposed over the last three decades, see Reference 4 for a review. The standard Tikhonov regularisations of zeroth, first or second order 3,5 are easy to implement thanks to the closedform solution of the quadratic inverse problem. Closely related to the Tikhonov regularisation is the generalised truncated singular value decomposition, in which small singular values of the forward operator are filtered out. 6 The inverse problem can also be interpreted from a Bayesian perspective. 7 In this case, the maximum a posteriori estimator, which matches the classical Tikhonov solution under the usual hypothesis of a standard Gaussian prior distribution, offers more flexibility in embodying prior knowledge in the inverse problem, for example, from a training data set. Nonsmooth regularisations such as total variation (TV) are a valid alternative to quadratic approaches in the presence of sharp gradients in the reconstructed data, 8 but lead to a significantly more difficult solution of the inverse problem. Approximated or smoothed versions of TV are therefore popular. 4 More recently, physiology-based and spatio-temporal regularisation approaches have also been considered. 9,10 The potential-based forward problem is particularly attractive when the torso is assumed as a homogeneous electric conductor. Then, the forward problem can be conveniently recast into an integral formulation involving only the boundaries of the torso, that is, the epicardium and the chest, with no need of solving the problem in the full threedimensional (3D) domain. 11,12 The numerical treatment of the boundary formulation is however less practical than a standard finite element approach in 3D, as it requires special care in the treatment of singular boundary integrals on piecewise smooth surfaces, like triangulated surfaces obtained from the segmentation of cardiac images. Acquired with given clinical constraints, imaging data are typically of limited resolution and noisy. Therefore, the segmentation of selected cavities is challenging and certainly subject to uncertainty. In the case of the heart, segmentation is made even more difficult by the movement of the organ. As a consequence, the resulting segmentation is subject to time-dependent shape uncertainty. This uncertainty in the anatomy propagates through the solution of the inverse problem, resulting in a reconstruction also affected by uncertainty.
In spite of its acknowledged importance, 13 uncertainty quantification in the context of cardiac modelling has emerged only very recently, see for example References 14-17. In electrocardiography, the most relevant uncertainty to account for is in the body surface electric recordings. Within a Bayesian framework, the inverse problem maps a prior distribution of the pericardial potentials into a posterior distribution, which also accounts for noise in the input data through the likelihood. 7 Model uncertainty has been considered as well. The electric conductivity of the torso is highly heterogeneous, for example, lungs, blood masses, interstitial tissue, muscles and bones significantly differ in terms of conduction, and uncertain, which may influence the reconstruction as well. 18,19 In the present context, shape uncertainty has however received very limited attention thus far, despite preliminary studies showed a non-negligible impact on the inverse reconstruction. 20,21 From the mathematical perspective, sophisticated tools and theory for the treatment of shape uncertainties are already available. Besides the fictitious domain approach considered in Reference 22, one typically distinguishes two approaches to deal with shape uncertainties: the perturbation method, see Reference 23 and the references therein, which is suitable to treat small perturbations of the nominal shape and the domain mapping method, see References 24-27. Here, we focus on the more flexible domain mapping approach. Then, the computation of quantities of interest, such as expectation and variance of the potential, gives rise to high dimensional quadrature problems for the random parameters. Given sufficient parametric regularity, sophisticated sparse quadrature and quasi-Monte Carlo methods, see for example References 28-32, can be applied. If such regularity is not present, one has to resort to the only slowly converging, Monte Carlo method.
This work aims at evaluating the impact of space-time uncertainty on the forward and inverse problems of electrocardiography. To model the shape uncertainty, we consider a space-time reference geometry and a random deformation field, that we represent by its Karhunen-Loève expansion. The covariance of the random deformation field depends on both space and time and accounts for the periodicity of the motion of the heart. Furthermore, we formulate the forward problem as a boundary integral equation by assuming a constant torso conductivity. We refer to Reference 12 for all the details of this approach and in particular for the applied discretisation. The boundary integral formulation is particularly suitable for our purpose, because both input data and observations are confined to a surface, these are the pericardium and the chest, respectively. In order to accelerate the computation of quantities of interest, we discretise the spatial problem with a rapidly converging collocation method, achieving high accuracy already with a relatively coarse mesh, while the deformation field is approximated by an efficient low-rank technique.
In order to assess the effect of the shape uncertainty, we estimate the expectation and the variance of the chest potential (forward problem) and the pericardial potential (inverse problem) using the anisotropic sparse quadrature method from Reference 28. This approach is validated against the quasi-Monte Carlo method based on Halton points, see Reference 32. The inverse problem, known to be severely ill-posed, needs an adequate regularisation. We employ classic regularisations proposed specifically for our problem of interest, such as zero order Tikhonov, the first-order Tikhonov and TV. 4 The implication of the choice of the regularisation on the convergence rate of the sparse quadrature is also analysed. Finally, we explore an H 1=2 regularisation, see References 33, a natural option stemming from the boundary integral formulation of the problem itself.
The paper is organised as follows. In section 2, we introduce the mathematical description of the problem at hand. Section 3 provides the corresponding boundary integral formulation and the discretisation of the involved boundary integral operators. The different regularisation methods for the inverse problem are presented in section 4. Section 5, dedicated to modelling shape uncertainty, is the cornerstone of the present work. Finally, section 6 is devoted to the numerical assessment of a two-dimensional, Cine Magnetic Resonance Imaging (MRI)-derived geometry. Herein, we validate the sparse quadrature and we quantify the effect of the shape uncertainty on the forward and inverse problems. Þdenote a complete and separable probability space, where Ω is a sample space, ℱ ⊆ 2 Ω is a σ-field and P : ℱ ! 0, 1 ½ is a probability measure. In what follows, we write Σ t,ω ð Þ with t ∈ 0, T ½ Þ and ω ∈ Ω to indicate the dependence of the pericardium on time and on the random parameter. Obviously, the time-dependence and the shape uncertainty of the heart also affect the torso, which we denote by D t, ω ð Þ, while we assume the chest Γ to be fixed over time and not being subject to uncertainty.
In the forward problem, given the potential u t ð Þ : Σ t, ω ð Þ!R at the pericardium, we wish to compute the potential on the chest Γ. Assuming that the torso is a homogeneous volume conductor and that no current sources are present, the extracellular potential y t, ω ð Þ: D t, ω ð Þ! in the whole torso satisfies the mixed boundary value problem We remark that in the model the electric conductivity is unit valued, with no loss of generality since its value, being constant, does not influence the solution. Moreover, we assume here that the potential u x, t ð Þ is given in spatial coordinates. A possibility to guarantee that u x, t ð Þ is well-defined for each realisation of the random parameter is to assume that it is given with respect to the hold-all domain which contains every possible space-time tube.

| The inverse problem of electrocardiography
Given the data y d t ð Þ : Γ ! R on the chest, the inverse problem corresponding to (1) is to find the potential u t, ω ð Þ: Σ t, ω ð Þ!R at the pericardium that satisfies This gives rise to the solution operator It is possible to show, see Reference 3, that A t,ω ð Þ is injective and continuous for t ∈ 0, T ½ and almost every ω ∈ Ω.
The minimum u t,ω ð Þ, however, is not stable with respect to perturbations. A common practice is to introduce a regularisation into the problem to recover stability, as follows: Herein, λ > 0 is the regularisation parameter and ℛ is the regularisation functional, for which several choices are possible. The optimum can be explicitly computed according to Later on we shall discuss different options for the regularisation.

| Boundary integral operators
In the following exposition, for the sake of an easier notation, we consider t and ω to be fixed. For a comprehensive exposition of the boundary integral approach, we refer the reader of Reference 12. For a better distinction of quantities that are defined with respect to the boundaries and those which are given on the domain, we introduce for Φ ∈ Σ, Γ, ∂D f gthe trace operators Then, we may write the forward problem according to For x ∈ D, the potential y x ð Þ is given by the representation formula where denotes the fundamental solution for the Laplacian. Introducing the single layer operator and taking the Dirichlet trace, Equation (2) yields the Dirichlet-to-Neumann map Next, introducing for Φ,Ψ ∈ Σ, Γ f g, the restricted operators we can split up Equation (4) into the system Rearranging this system in order to move the unknowns γ int 1,Σ y and γ int 0,Γ y to the left side and the data γ int 0,Σ y ¼ u and γ int 1,Γ y ¼ 0 to the right side, we arrive at the system boundary integral equations As V ΣΣ is an elliptic operator, considering the Schur complement yields the operator equation for the desired potential γ int 0,Γ y on the chest. In particular, the solution operator A : while the Poincaré-Steklov operator ℬ :

| Numerical discretisation
Our goal is to discretise the boundary integral equations derived previously. In order to not having to deal with nonconstant coefficients, we shall solve the equations in the spatial coordinate frame, that is, for each tuple t, ω ð Þ, assemble the boundary integral operators on Σ t, ω ð Þ and Γ. To simplify the notation, in the following quantities we omit the dependency of Σ on t and ω. As the focus of this work is on uncertainty quantification, rather than on the solution of the spatial problem, we restrict ourselves to the simplified 2D anatomy of Figure 1. Hence, from (3), we have Given that the domain D t, ω ð Þ exhibits C 2 -boundaries, which are described by the parameterisations γ Σ : 0, 1 ½ Þ! Σ t, ω ð Þ and γ Γ : 0, 1 ½ Þ!Γ, we can employ the collocation method from Reference 12. This method is based on the trapezoidal rule and an appropriate desingularisation of the kernel functions. In view of the Euler-Maclaurin formula, it converges quadratically for C 2 -boundaries and exponentially for smooth boundaries. More precisely, given that the solution to (5) satisfies y ∈ C k ∂D ð Þ, there holds for some C > 0, where y n is obtained from the collocation method with n ¼ 2j points for some j ∈ ℕ. We refer to Reference 12 for all the details. For the sake of completeness, we briefly recall the collocation method in Appendix A.

| SOLUTION OF THE INVERSE PROBLEM
In the inverse problem, the goal is to reconstruct the potential u at the pericardium Σ t, ω ð Þ given noisy data measurements y d ∈ H 1=2 Γ ð Þ on the chest. In the following exposition, for the sake of an easier notation, we indicate the dependency of Σ on the tuple t, ω ð Þ only when specifying function spaces. In particular, we consider the solution operator From (8), we deduce that there exist a matrix A ∈ R n Γ Ân Σ and a matrix B ∈ R n Σ Ân Σ such that ρ 0,Γ ¼ Au and ρ 1,Σ ¼ Bu: In this section, we consider inverse problem solutions of the form where λ > 0 is the regularisation parameter and ℛ is the regularisation functional. We discretise (9) employing, for The accuracy of this approximation is consistent with the one obtained by the collocation method, as can easily be inferred from the Euler-Maclaurin formula. We arrive at the discrete formulation where R is the discrete regularisation term. The optimisation of the objective function leads to the equation for the minimum u. In the following subsections, we present four strategies to solve the inverse problem based on different regularisations.

| Zero order Tikhonov regularisation
The zero order Tikhonov regularisation penalises the L 2 -norm of the Dirichlet data at Σ t,ω ð Þ, see Reference 34, and the regularisation functional in (9) reads The discrete regularisation term in (10) is thus and in (11) we have

| First-order Tikhonov regularisation
The first-order Tikhonov regularisation penalises the L 2 -norm of the Neumann data at Σ t, ω ð Þ, see Reference 34, and the regularisation functional in (9) is Then the discrete regularisation term in (10) is and in (11) we have

| H 1=2 regularisation
When penalising the H 1=2 -norm of the Dirichlet data at Σ t, ω ð Þ, the regularisation functional in (9) is given by The resulting discrete regularisation term in (10) is and its Fréchet derivative in (11) is

| Total variation regularisation
In the TV regularisation, we penalise the L 1 -norm of the Neumann data at Σ t,ω ð Þ, that is, the regularisation functional in (9) is As the L 1 -norm is non-differentiable, it is common to employ the approximation where β > 0 is a small constant (typically β ¼ 10 À5 ). This introduces a non-linearity in the optimality condition. This difficulty is handled in Reference 4 by using the zero order Tikhonov solution u 0 in the non-linear term. This approach is equivalent to set (9). The discretisation requires the definition of the diagonal matrix W β u 0 ð Þ∈ R n Σ Ân Σ with entries for i ¼ 0, …, n Σ À 1, where u 0 is the discrete zero order Tikhonov solution. The discrete regularisation term in (10) then reads

| Modelling time-dependent shape uncertainty
In order to assess the shape uncertainty of the pericardium Σ t, ω ð Þ, we assume the existence of a reference domain D ref & R d , with boundaries Σ ref and Γ, and a random deformation field Herein, the expectation is given in terms of the Bochner integral we can interpret χ as a random deformation field χ b z, ω ð Þ acting on the space-time cylinder. A similar model has originally been introduced in Reference 25 for the stationary case. Assuming it is possible to expand χ b z,ω ð Þ in a Karhunen-Loève expansion, cf., Reference 35 according to Herein, λ k ,χ k ð Þare the eigen-pairs of the integral operator defined by the matrix valued covariance function while Y k f g is a family of uncorrelated random variables with normalised variance. To avoid distortion of the space-time tubes χ Q ref ,ω ð Þ, ω ∈ Ω, and hence to guarantee well-posedness of the boundary value problem at hand, cf. (1), we impose the uniformity condition for some constant C uni > 0 uniformly in ω ∈ Ω. 1 As a consequence, the random variables' Y k f g ranges are bounded. Therefore, without loss of generality we may assume Y k : Ω ! À1, 1 ½ , k ¼ 1, 2, …. In particular, we obtain Replacing the random variables by their image À1, 1 ½ ℕ , we arrive at the parametrised Karhunen-Loève expansion For the numerical computation of the Karhunen-Loève expansion, it is sufficient to know the expectation and the covariance at Σ ref Â 0, T ½ , see Reference 36. In our concrete case, it is even sufficient to only compute the Karhunen-Loève expansion in the space-time collocation points. To this end, we employ the pivoted Cholesky decomposition, see References 25,37,38. These results in a finite rank Karhunen-Loève expansion where b z i are the space-time collocation points on Σ ref Â 0, T ½ and K is the dimension of the random parameter space.

| Computation of quantities of interest
We make the common assumption that the random variables are even independent and uniformly distributed. Then, the push-forward measure is given by the product density Hence, we can now express quantities of interest, such as the moments of the potential on the chest, by means of an integral with respect to À1, 1 ½ K . It holds for the moments that Especially, the expectation and the variance are given by In practice, the moments need to be approximated by a quadrature rule, that is, where ξ i ∈ À1, 1 ½ K for i ¼ 1, …, N are the quadrature points and w i ∈ R are the corresponding weights. Due to the spacetime modelling, the random parameter space is very high-dimensional, and appropriate quadrature rules need to be employed. In this article, we employ the anisotropic sparse quadrature from Reference 28, which is based on the Gauss-Legendre quadrature rule. For the sparse quadrature to converge, the quantity of interest needs to be regular with respect to the random parameter ξ ∈ À1, 1 ½ K . For the forward problem, such results are available, we refer to Reference 25. For the Dirichlet control problem with an affine random diffusion coefficient, such results have also been derived in Reference 39.

| NUMERICAL EXPERIMENTS
In this section, we present numerical experiments to assess the validity of the approach. To obtain a space-time reference torso anatomy, we have segmented cardiac magnetic resonance (CMR) imaging data previously acquired. 40 The temporal image stack was obtained from a cine ECG-triggered segmented steady-state free precession sequence in midventricular short-axis orientation. Slice thickness and voxel resolution were 8 and 1.36 mm, respectively. The 25-phase temporal stack covered the whole heartbeat, with an RR interval on the surface ECG of T ¼ 690 ms. Images were reordered such that the first image at t ¼ 0 ms corresponded to the diastole, defined by the maximal left ventricular cavity volume. The systole, defined by the minimum cavity volume, occurred at t ¼ 270 ms. The segmentation was performed by manual contour tracing of the epicardium for each image. In order to end up with a smooth computational domain, we performed a least-squares fit of the contours using a truncated Fourier series with a threshold of 10 À3 in the relative root-mean-square error. Finally, we interpolated the extracted shapes to get a pericardial representation at 50 time instants. The chest was also previously segmented from ultra-fast gradient-echo "VIBE" images in axial, coronal, and sagittal orientations to produce a smooth three-dimensional closed surface modelled in Blender.2 As shown in Figure 1, the contour of the chest was eventually obtained by intersecting the chest surface with the orientation plane of the pericardium.
In Figure 2, we summarised the input data for the forward and inverse problem and the geometry of the pericardium, superimposed with the CMR images. To model the uncertainty, we selected the resulting reference shape of the pericardium Σ ref t ð Þ, for t ∈ 0, T ½ Þ, as the space-time mean of the random deformation field χ . The covariance of χ evaluated at the points Þ was the product of a matrix-valued kernel in space and a scalar kernel in time. For the spatial component of the covariance, we considered the Matérn kernel with parameters ρ ¼ 50, σ 2 ¼ 4=3 and different values of the smoothness index ν. Note that in the definition of k M ν , the function Γ Á ð Þ is the gamma function and K ν Á ð Þ is the modified Bessel function of the second kind. Since we are interested in modelling the uncertainty due to the noise and the limited resolution of CMR, the distance x 0 used in the definition of the kernel, was measured with the Euclidean norm in R 2 . In contrast, for the distance in time between t and t 0 , we take into account the periodicity of the motion of the pericardium. To do so, we first scale the interval 0,T ½ Þ to 0, 2π ½ Þwith the mapping F I G U R E 2 Geometry and input data. First row: cardiac magnetic resonance images with superimposed segmented chest (blue) and time-dependent pericardium (red). The shaded region around the pericardium represents the shape confidence interval, obtained as  Σ ½ t ð ÞAE1:96 Á StdDev Σ ½ t ð Þ. The yellow dots correspond to γ Γ 0 ð Þ and γ Σref t ð Þ 0 ð Þ. Second row: forward data u γ Σ s ð Þ, t ð Þ . Third row: inverse data y d γ Γ s ð Þ, t ð Þ and then we compute the geodesic distance θ t,t 0 ¼ arccos cosðη t ð Þ Àηðt 0 ÞÞ : Finally, to model an additional correlation between neighbouring time slices, we employ the sine power kernel see Reference 41. Employing a tensor product construction, we end up with the covariance function

5:
Note that we assume here that there is no correlation between the two spatial components of the deformation field. By construction, the joint covariance kernel is positive-definite on the space R 2 Â  1 , where  1 is the unit circle.
In the forward problem, the pericardial potential u x,t ð Þ was defined analytically as a T-periodic function in the variable t. For convenience, we set u x,t ð Þ¼u γ Σ t,ω ð Þ s ð Þ, t , being s ∈ 0, 1 ½ Þ the normalised curvilinear coordinate. We simulated a left bundle branch block, that is the extracellular potential consisted in a propagation from free wall of the right ventricle, at s ¼ 0, towards the free wall of the left ventricle, at s ¼ 0:5. The propagation took 150 ms, consistent with a long QRS complex. The specific analytical form of u x, t ð Þ is given in Appendix B. To generate the input data for the inverse problem, we computed the solution of the forward problem on the spacetime reference geometry and eventually added Gaussian noise with zero mean and variance 10 À8 , corresponding to a signal-to-noise ratio of approximately 46 dB.
Concerning the space discretisation, we considered n Σ ¼ n Γ ¼ 500 collocation points. This resulted in a truncated Karhunen-Loève expansion with K ¼ 648 terms, obtained with a tolerance of 10 À4 . That is, the parametric dimension of the UQ problem was 648, thus high-dimensional. We remark that, since the boundaries are represented by trigonometric polynomials and all data are smooth functions, the collocation method converges exponentially. Resulting in spatial approximation errors for the chosen number of boundary points, which are already of the order of the machine precision.
Since we wish to employ the sparse quadrature to estimate our quantities of interest, we first numerically tested its convergence by comparing it to the quasi-Monte Carlo method based on the Halton set, see Reference 32. We tested the F I G U R E 3 Convergence plot of the first (ℳ 1 ) and second (ℳ 2 ) moment of the forward solution convergence of the first and second moment of the forward and inverse solution at t ¼ 189 ms. The time-independent problem yielded a truncated Karhunen-Loève expansion of K ¼ 101 terms. To validate the applicability of the sparse quadrature, we consider the approximation error of the moments by a comparison to the quasi-Monte Carlo method based on Halton points. The accuracy of approximately 5 Á 10 À5 was achieved by using 17'799 quadrature points within the sparse quadrature. Figure 3 shows the convergence plot of the forward solution. We experimentally observed a convergence rate of 0.75. The obtained result indeed corroborates that the forward problem is smooth in the parametric space.
For the inverse problem, we tested the regularisations reported in section 4. The regularisation parameter was determined with the L-curve method, see Reference 34, on the space-time reference geometry. For each regularisation, we then chose the maximal regularisation parameter over all the time steps. This might have led to over-regularisation for some time steps, but it has the benefit of limiting the oscillations in the numerical solution of the inverse problem. The values of the chosen regularisation parameters λ are reported in Table 1.  Figure 4 shows the convergence plots of the inverse solutions. We can observe convergence towards the reference moments for the zero order Tikhonov, the first-order Tikhonov and the H 1=2 regularisations. Instead, the curves resulting from the TV regularisation flatten out, meaning that there is no convergence and suggesting that this regularisation is not sufficiently smooth with respect to the random parameter, and therefore inadequate for the sparse quadrature approach. In contrast, the other regularisations are to be suitable for the sparse quadrature approach. Next, we present the corresponding numerical results for the inverse problem, when employing the H 1=2 regularisation. We estimate the mean and the standard deviation of the forward and inverse solution. The computed quantities of interest for the forward problem are shown in Figure 5. We can observe that the expectation was very close to the chest potential computed on the space-time reference geometry. Moreover, the standard deviation is small and it mainly affects the regions with higher amplitude during the depolarisation phase. In particular, the standard deviation is higher close to γ Γ 0 ð Þ. The reason for this phenomenon is that points in the vicinity of γ Γ 0 ð Þ are close to the heart, see Figure 2. Therefore the effect of the shape uncertainty on the forward solution is limited.
The computed quantities of interest for the inverse problem are shown in Figure 6. We observe that the expectation deviates slightly from the pericardial potential and shows some oscillations in the depolarisation phase. This is due to the ill-posedness of the inverse problem and the effect of the regularisation. Indeed, the Tikhonov regularisation, especially at lower order, tends to introduce oscillations in the solution, if the regularisation parameter is not sufficiently large. On the other hand, a strong regularisation yields smoothed out gradients, thus less accurate inverse solution. The H 1=2 regularisation might therefore introduce oscillations as well. In the bottom row of Figure 6, the standard deviation is larger than in the forward problem and the shape uncertainty mainly affects regions with large gradients. In the depolarisation phase, the potential is steeper and, consequently, the standard deviation is larger. This fact may have a consequence on the quality of derived quantities such as the activation map, usually obtained from the point with largest negative deflection in the signals. In conclusion, the effect of the shape uncertainty is much more relevant for the inverse problem than for the forward problem.

| CONCLUSIONS
In this work, we have considered the forward and inverse problem of electrocardiography in the presence of space-time shape uncertainties. The high parametric dimensionality of the problem, along with the non-linear map between the input and output uncertainties, renders the problem challenging from the mathematical and the numerical perspective.
To address space-time shape uncertainties, we have suggested a model by random space-time deformation fields with a periodic-in-time covariance kernel. This approach resulted in a high-dimensional uncertainty quantification problem. To make this dimensionality feasible for numerical computations, we rely on the one hand on a low-rank representation of the random deformation field and efficient quadrature techniques on the other hand. The resulting spatial problem for each quadrature point in the parameter was addressed by a boundary integral formulation in combination with a highly accurate and fast collocation method. The numerical results indicate that shape uncertainties affect the solution much more strongly in the inverse problem than in the forward one, as a consequence of the ill-posedness. Moreover, it was observed that the regularisation of the inverse solution may affect the regularity of the solution with respect to the stochastic parameter, with possible limitations in the convergence order of the quadrature method. Especially, the TV regularisation showed poor performance in this respect, while the H 1=2 regularisation is suitable for this problem, both in terms of regularity and quality of the reconstruction.
Future research directions include the transition to 3D models and the generalisation to piecewise-constant conductivities in the torso. In this case, it is still possible to leverage on the boundary integral formulation. From the mathematical perspective, a general result on the parametric regularity of the inverse problem is still missing. Such a result would help to shed some light on the sub-optimal convergence rate observed in our experiments and the apparent lack of parametric regularity of the inverse problem with TV regularisation. Finally, from the clinical perspective, we shall further investigate how segmentation uncertainty may be modelled. As a matter of fact, the uncertainty likely results from multiple sources. Besides noise, CMR images are affected by breath holding, which may limit the volume of the heart and spatio-temporal mis-registration. The segmentation process may amplify the acquisition uncertainty, depending on the method. Importantly, the segmentation is sometimes a non-smooth operation, in the sense that small perturbation in the images could yield large and possibly topological changes in the contour. Therefore, special care will be needed in order to correctly model the uncertainty from clinical data.