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

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 H 1 / 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. ECG imaging technologies have been extensively validated in experimental, animal and, more recently, clinical settings, with promising results.

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. 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 , are easy to implement thanks to the closed‐form 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. The inverse problem can also be interpreted from a Bayesian perspective. 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. Non‐smooth regularisations such as total variation (TV) are a valid alternative to quadratic approaches in the presence of sharp gradients in the reconstructed data, but lead to a significantly more difficult solution of the inverse problem. Approximated or smoothed versions of TV are therefore popular. More recently, physiology‐based and spatio‐temporal regularisation approaches have also been considered. ,

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 three‐dimensional (3D) domain. , 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, uncertainty quantification in the context of cardiac modelling has emerged only very recently, see for example References 14, 15, 16, 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. 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. , 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. ,

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, 25, 26, 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, 29, 30, 31, 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. The implication of the choice of the regularisation on the convergence rate of the sparse quadrature is also analysed. Finally, we explore an 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.

PROBLEM FORMULATION

The forward problem of electrocardiography

Mathematically, we start from a smooth domain representing the torso, whose boundaries are the chest denoted by and the pericardium denoted by . In other words, we have and , see Figure 1 for a visualisation of the domain.

FIGURE 1

2D cross section of the torso reconstructed from an MRI dataset

The pericardium depends periodically on time and is subject to uncertainty. To model the time‐dependence, let be the time interval of interest, where is the duration of one heartbeat. To model the uncertainty, let denote a complete and separable probability space, where is a sample space, is a ‐field and is a probability measure. In what follows, we write with 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 , while we assume the chest to be fixed over time and not being subject to uncertainty.

In the forward problem, given the potential 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 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 is given in spatial coordinates. A possibility to guarantee that 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 on the chest, the inverse problem corresponding to (1) is to find the potential at the pericardium that satisfies

In the forward problem, if , then the solution satisfies and its trace on is in . This gives rise to the solution operator

It is possible to show, see Reference 3, that is injective and continuous for and almost every . In the inverse problem, therefore, if then there exists a unique minimum such that at . The minimum , 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, 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 FORMULATION

Boundary integral operators

In the following exposition, for the sake of an easier notation, we consider 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 the trace operators

Then, we may write the forward problem according to

For , the potential is given by the representation formula where denotes the fundamental solution for the Laplacian. Introducing the single layer operator the double layer operator and taking the Dirichlet trace, Equation (2) yields the Dirichlet‐to‐Neumann map

Next, introducing for , the restricted operators and we can split up Equation (4) into the system

Rearranging this system in order to move the unknowns and to the left side and the data and to the right side, we arrive at the system boundary integral equations

As is an elliptic operator, considering the Schur complement yields the operator equation for the desired potential on the chest. In particular, the solution operator is given by while the Poincaré–Steklov operator is given by

Numerical discretisation

Our goal is to discretise the boundary integral equations derived previously. In order to not having to deal with non‐constant coefficients, we shall solve the equations in the spatial coordinate frame, that is, for each tuple , assemble the boundary integral operators on and . To simplify the notation, in the following quantities we omit the dependency of on 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 and

Given that the domain exhibits ‐boundaries, which are described by the parameterisations and , 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 ‐boundaries and exponentially for smooth boundaries. More precisely, given that the solution to (5) satisfies , there holds for some , where is obtained from the collocation method with points for some . We refer to Reference 12 for all the details. For the sake of completeness, we briefly recall the collocation method in Appendix A.

For , we consider the ‐points discretisations for . Moreover, we denote by the identity matrix. The discretisation results in the block linear system cp. (5). Herein

Furthermore, we set and is the ‐dimensional vector of zeros.

SOLUTION OF THE INVERSE PROBLEM

In the inverse problem, the goal is to reconstruct the potential at the pericardium given noisy data measurements on the chest. In the following exposition, for the sake of an easier notation, we indicate the dependency of on the tuple only when specifying function spaces. In particular, we consider the solution operator , cp. (6), and the Poincaré–Steklov operator , cp. (7). In the discrete setting, we define the vectors and

From (8), we deduce that there exist a matrix and a matrix such that

In this section, we consider inverse problem solutions of the form where is the regularisation parameter and is the regularisation functional. We discretise (9) employing, for , the trapezoidal rule where is the diagonal mass matrix with entries 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 is the discrete regularisation term. The optimisation of the objective function leads to the equation for the minimum . 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 ‐norm of the Dirichlet data at , 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 ‐norm of the Neumann data at , see Reference 34, and the regularisation functional in (9) is

Then the discrete regularisation term in (10) is and in (11) we have

regularisation

When penalising the ‐norm of the Dirichlet data at , 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 ‐norm of the Neumann data at , that is, the regularisation functional in (9) is

As the ‐norm is non‐differentiable, it is common to employ the approximation where is a small constant (typically ). This introduces a non‐linearity in the optimality condition. This difficulty is handled in Reference 4 by using the zero order Tikhonov solution in the non‐linear term. This approach is equivalent to set in (9). The discretisation requires the definition of the diagonal matrix with entries for , where is the discrete zero order Tikhonov solution. The discrete regularisation term in (10) then reads while in (11) it results in

UNCERTAINTY QUANTIFICATION

Modelling time‐dependent shape uncertainty

In order to assess the shape uncertainty of the pericardium , we assume the existence of a reference domain , with boundaries and , and a random deformation field such that where

Herein, the expectation is given in terms of the Bochner integral

Introducing the space–time cylinder with coordinates we can interpret as a random deformation field 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 in a Karhunen–Loève expansion, cf., Reference 35 according to

Herein, are the eigen‐pairs of the integral operator defined by the matrix valued covariance function while is a family of uncorrelated random variables with normalised variance.

To avoid distortion of the space–time tubes , , and hence to guarantee well‐posedness of the boundary value problem at hand, cf. (1), we impose the uniformity condition for some constant uniformly in . As a consequence, the random variables' ranges are bounded. Therefore, without loss of generality we may assume , . In particular, we obtain

Replacing the random variables by their image , 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 , 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 are the space–time collocation points on and 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 . 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 for are the quadrature points and are the corresponding weights. Due to the space–time 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 . 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. The temporal image stack was obtained from a cine ECG‐triggered segmented steady‐state free precession sequence in mid‐ventricular 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 ms. Images were reordered such that the first image at ms corresponded to the diastole, defined by the maximal left ventricular cavity volume. The systole, defined by the minimum cavity volume, occurred at 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 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 , for , as the space–time mean of the random deformation field . The covariance of evaluated at the points with and with 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 , and different values of the smoothness index . Note that in the definition of , the function is the gamma function and 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 used in the definition of the kernel, was measured with the Euclidean norm in . In contrast, for the distance in time between and , we take into account the periodicity of the motion of the pericardium. To do so, we first scale the interval to with the mapping and then we compute the geodesic distance

FIGURE 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 . The yellow dots correspond to and . Second row: forward data . Third row: inverse data

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

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 , where is the unit circle.

In the forward problem, the pericardial potential was defined analytically as a ‐periodic function in the variable . For convenience, we set , being 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 , towards the free wall of the left ventricle, at . The propagation took 150 ms, consistent with a long QRS complex. The specific analytical form of is given in Appendix B.

To generate the input data for the inverse problem, we computed the solution of the forward problem on the space–time reference geometry and eventually added Gaussian noise with zero mean and variance , corresponding to a signal‐to‐noise ratio of approximately 46 dB.

Concerning the space discretisation, we considered collocation points. This resulted in a truncated Karhunen–Loève expansion with terms, obtained with a tolerance of . That is, the parametric dimension of the UQ problem was , 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 convergence of the first and second moment of the forward and inverse solution at  ms. The time‐independent problem yielded a truncated Karhunen–Loève expansion of 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 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.

FIGURE 3

Convergence plot of the first () and second () moment of the forward solution

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.

TABLE 1

Regularisation parameter choice for the zero‐order Tikhonov, first‐order Tikhonov, and total variation regularisations

	Zero‐order Tikhonov	First‐order Tikhonov	H^1/2	Total variation
λ	10^−6	10^−3	10^−5	10^−5

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 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.

FIGURE 4

Convergence plot of the first () and second () moment of the inverse solution with (A) zero‐order Tikhonov, (B) first‐order Tikhonov, (C) and (D) total variation regularisation

Next, we present the corresponding numerical results for the inverse problem, when employing the 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 . The reason for this phenomenon is that points in the vicinity of are close to the heart, see Figure 2. Therefore the effect of the shape uncertainty on the forward solution is limited.

FIGURE 5

Solution of the forward problem. First row: expectation and confidence interval (blue), and solution with reference geometry (dashed black). Second row: contour plots in space–time of the reference solution, expected solution and standard deviation

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 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.

FIGURE 6

Solution of the inverse problem. First row: expectation and confidence interval (red), and reference solution (dashed black). Second row: contour plots in space–time of the reference solution, expected solution and standard deviation

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 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.