# A Volume Source Method for Solving ECGI Inverse Problem

Mohamadou Malal Diallo, Yves Coudière, Rémi Dubois

► To cite this version:

Mohamadou Malal Diallo, Yves Coudière, Rémi Dubois. A Volume Source Method for Solving ECGI Inverse Problem. FIMH 2021 - 11th International Conference on Functional Imaging and Modeling of the Heart, Jun 2021, Stanford / Virtual, United States. pp.551-560, 10.1007/978-3-030-78710-3_53 .
hal-03507452

HAL Id: hal-03507452

https://hal.science/hal-03507452

Submitted on 3 Jan 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# A volume source method for solving ECGI inverse problem

Mohamadou Malal Diallo<sup>1,2,3</sup>, Yves Coudière<sup>1,2,3</sup>, and Rémi Dubois<sup>1,4,5</sup>

<sup>1</sup> IHU-Liryc, Fondation Bordeaux Université, Pessac, France

<sup>2</sup> CARMEN Research Team, Inria Bordeaux – Sud-Ouest, Talence, France

<sup>3</sup> Univ. Bordeaux, IMB, UMR 5251, CNRS, Talence, France

<sup>4</sup> Univ. Bordeaux, CRCTB de Bordeaux, U1045, Bordeaux, France

<sup>5</sup> INSERM, CRCTB, U1045, Bordeaux, France

mohamadou-malal.diallo@inria.fr

**Abstract.** Electrocardiographic Imaging (ECGI) is a non-invasive procedure that allows to reconstruct the electrical activity of the heart from body surface poten- tial map (BSPM). In this paper, we present a volume model to solve the elec- trocardiography inverse problem capable to take into account structural informa- tions obtained by imaging techniques. The direct problem maps a volume current in the cardiac muscle (ventricles) to the body surface electrical measures. The model is based on coupling bidomain heart model with torso conduction. The corresponding inverse problem is solved with the Tikhonov regularization. Sim- ulated database are used for the evaluation of this method and we compared them to standard method of fundamental solutions (MFS). The sensitivity to noise is also assessed. The correlation coefficients (CC) and the relative error (RE) of activation times were computed. Results show that the CC (respectively RE) me- dian is respectively 0.75 for the volume model and 0.4 for the MFS (respectively 0.31 vs 0.35) on the epicardium. On the endocardium, the CC and the RE median are 0.65 and 0.33 for the volume method. In conclusion, the volume method per- forms better than the method of fundamental solutions (MFS) for any noise level, and reconstruct in addition endocardial information.

**Keywords:** Inverse problem · Scar · ECGI · Electrocardiography.

## 1 Introduction

Electrocardiographic imaging (ECGI) is a non-invasive technique that is used to re-construct the electrical activity of the heart from body surface electrical potential maps (BSPM), and the geometry of the heart and torso.

The most common approach to compute this reconstruction follows the electrostatic theory of perfect volume conductors. It assumes that the electrical field outside of the heart (or just the ventricles in our case) is perfectly defined by the outward current flux or potential distribution on its surface, considered as enclosing a finite volume. In such an approach, the distribution of charges in the enclosed volume is not relevant. Hence the Laplace equation is set on the volume between the epicardial and the body surfaces, and defines a direct mapping from epicardial potentials to body surface potentials. The
inverse problem consists in finding the epicardial potentials from the body surface potentials, assuming the additional zero flux condition on the body surface. This is a famous Cauchy problem, notably ill-posed [1]. It is commonly solved as an optimization problem with a Tikhonov regularization. The method of fundamental solutions (MFS) with Tikhonov regularization is commonly used in this case [8]. In clinical care, structural images of the patient are often available, that may show details about the electrical properties inside the cardiac volume. A major question is therefore how to integrate this information in order to drive the inverse problem. The classical approach is limited because it cannot take into account scars inside the heart volume, and because it looks for the electrical field on the epicardial surface.

The homogenized bidomain equations describe the electrical behavior of the heart within the torso in a different manner, since they consider two superimposed electrical fields, the intracellular one $u_i$ supposed to be defined in the myocardium only, with a no flux boundary condition, and the extracellular one $u_e$, that prolong as the extracardiac one outside the heart. The homogenized equations are strongly coupled, since both electrical fields are found to be defined by an anisotropic Laplace equation with a right hand side equal to $\xi (C_m \partial_t (u_i - u_e) + I_{\text{ion}})$, that is the transmembrane current per unit of surface rescaled by the ratios of surface of membrane per unit volume of tissue. This right-hand side is a current per unit of volume in A.m$^{-3}$.

In this paper, we consider the torso as a heterogeneous volume conductor for the electrical potential $u$ equal to the extracellular potential field in the myocardium, and to the extracardiac one elsewhere. This volume conductor includes the intracavitary blood, the heart and the remaining torso volume. The electrical potential $u$ is controlled by the data of the current per unit of volume defined above. We recall that a similar model is used in other studies like [5, 7], where the authors looked to reconstruct the transmembrane potential (TMP) source in the specific case of ischemia.

Our direct problem maps the current per unit of volume, called *volume current source*, in the heart to the BSPM. We can hence try to reconstruct the volume current source, and we are able to take into account different electrical conductivity values in each of the regions, and in particular in scars. We use the approach based on transfert matrix for solving the inverse problem, and it is computed by using a standard finite element method (FEM).

The problem is generally ill-posed and underdetermined, hence we apply the Tikhonov regularization technique. The regularization applies to the volume current source, while it applies on virtual charges spread around the torso volume for the MFS. We will refer to this method as the volume method (VM). We wanted to test the method on public datasets, and therefore used the ECGI consortium database EDGAR [6]. Among EDGAR datasets, only one contains enough data to test our method (KIT-20-PVC-Simulation-1906-10-30-EP-EndoEpi). The dataset is computed on realistic human-based anatomical model. The activation times (AT) recovered by the standard MFS and VM method were compared to the reference AT obtained from the model. The noise sensitivity of the model was also assessed.
# 2 Methods

## 2.1 Mathematical model

The mathematical model derives from the standard bidomain equations, and has the form of a Laplace equation for the potential field $u(t,x)$ in the blood, heart, and torso domains, respectively denoted by $\Omega_B$, $\Omega_H$, and $\Omega_T$ (figure 1). It reads

$$-\nabla \cdot (\sigma(x)\nabla u(t,x)) = F(t,x) \quad \text{in } \Omega, \qquad (1)$$

(where $\Omega = \Omega_B \cup \Omega_H \cup \Omega_T$) with the no flux boundary condition $-\sigma\nabla u \cdot n = 0$ on the torso surface $\partial\Omega$. In this model, the electrical conductivity is the function $\sigma(x)$, piecewise constant:

$$\sigma(x) = \begin{cases} \sigma_H & \text{if } x \in \Omega_H \setminus \Omega_{\text{scar}}, \\ \rho \sigma_H & \text{if } x \in \Omega_{\text{scar}}, \\ \sigma_B & \text{if } x \in \Omega_B, \\ \sigma_T & \text{if } x \in \Omega_T, \end{cases} \qquad (2)$$

where the factor $\rho \le 1$ may be used to decrease the conductivity in the scar volume $\Omega_{\text{scar}} \subset \Omega_H$ (possibly $\Omega_{\text{scar}} = \emptyset$).

The source term $F(t,x)$ accounts for the total ionic and diffusion currents, and has the form

$$F(t,x) = \begin{cases} f(t,x) & \text{if } x \in \Omega_H, \\ 0 & \text{if } x \in \Omega_B \cup \Omega_T. \end{cases} \qquad (3)$$

Thus, the geometry and the conductivity coefficients are parameters of our model. The

<figure><img src="image_1.png" /><figcaption>Fig. 1. Sketch of the model with a scar.</figcaption></figure>

equation derives from the bidomain one, and the potential field $u(t,x)$ is interpreted as the extracellular in the heart $\Omega_H$ and the extracardiac potential in the rest of the domain. The body surface potential map is the trace of $u(t,x)$ on the boundary of $\Omega$, denoted by $u_T := u|_{\partial\Omega}$. At each time instant, the forward problem maps the volume current source $f(t,x)$ to the BSPM $u_T(t,x)$. The inverse problem consists in finding an inverse mapping from BSPM data $d$ to a volume current source $f$ in the heart, which direct image best matches $d$.

In the present work, we put the factor $\rho = 1$ because there is no scar area in dataset.
## 2.2 Numerical resolution

Equation (1) was discretized by the P1 Lagrange finite element method (FEM) using the FreeFem++ software [4]. The matrix of the direct problem, called the transfer matrix, is an $N_T \times N_H$ matrix denoted by $M$, such that $u_T = Mf$, where $N_T$ is the number of electrode nodes on the body surface $\partial\Omega$, and $N_H$ the number of mesh nodes in the heart domain $\Omega_H$. The problem is generally ill-posed and underdetermined, and it need to be regularized for instance by a Tikhonov term. Given data $d \in \mathbb{R}^{N_T}$ on electrode torso nodes, we reconstructed the volume current source $f \in \mathbb{R}^{N_H}$ as the optimum

$$f = \operatorname*{argmin}_{g \in \mathbb{R}^{N_H}} \left\{ \| Mg - d \|_{l_2(\mathbb{R}^{N_T})}^2 + \lambda^2 \| g \|_{l_2(\mathbb{R}^{N_H})}^2 \right\}, \quad (4)$$

where $\lambda$ is the regularization parameter. The solution of problem (4) is finally computed by solving the normal equation

$$(M^T M + \lambda^2 I) f = M^T d, \quad (5)$$

where $I$ is the identity matrix. The CRESO method was used to choose the regularisation parameter $\lambda$ and the median of all time step regularization parameter was chosen for global parameter.

Given a BSPM time sequence $d = (d_i^n)_{i \in \{0, \dots, N_T\}}$ with $n \in \{0, \dots, T\}$, we finally have reconstructed a volume current source time sequence $f = (f_j^n)_{j \in \{0, \dots, N_H\}}$, and are able to reconstructed an extracellular time sequence $(u_j^n)_{j \in \{0, \dots, N_H\}}$ by solving equation (1). The reference MFS gives acces a sequence of potentials on the epicardium (with nodes $N_e$) $(w_k^n)_{k \in \{0, \dots, N_e\}}$. For each sequence, we define the activation times (AT) as the time of maximum negative derivation along time. It reads as follows :

$$AT_i = \operatorname*{argmin}_{n \in \{0, \dots, T\}} (s_i^{n+1} - s_i^n),$$

where $s_i^n$ is the signal time sequence at each point $P_i$ on the surface (epicardium or endocardium) at time $n$. We used also the method proposed by Duchateau et al [3] to smooth the AT maps. It consists to estimate the delays in activation for neighboring signal locations.

## 2.3 Simulated data and evaluation

Datasets were provided by the ECGI consortium database EDGAR, and were simulated in Karlsruhe Institute of Technology (KIT) [6]. The simulations pacing were performed on a voxel-based grid using a cellular automaton, then the BSPM of these simulated beats were interpolated on tetrahedral mesh and forward-calculated using a finite element method and a bidomain model. The extracellular potentials were extracted from the tetrahedra mesh at 163 electrode positions and at the nodes of the endo- and epicardial surface mesh. The heart model consisted of left and right ventricles and it was based, with thorax model, on real human anatomy MRI data. The data include 8 sets obtained from simulations of single pacing sites distributed as follows : center of septum (SEPTCENTER), left ventricle apex (LVAPEX), left ventricle lateral (LVLAT), left
ventricle anterior (LVANT), left ventricle lateral epicardium (LVLATEPI), left ventricle lateral endocardium (LVLATENDO), right ventricle posterior (RVPOST) and right ventricle anterior (RVANT).

Reference activation times (AT) were computed directly from the value of simulated extracellular potential. Using the simulated BSPM, the volume current source were computed with the VM. Afterwards AT were calculated using the method proposed by Duchateau [3] directly for the volume current source. Using the same data, the extracellular epicardial potentials were computed with the classical MFS, and AT calculated.

To evaluate the methods, we compute the correlation coefficients (CC) and the relative error (RE) on the activation times (AT) defined as follows :

$$CC = \frac{\sum_{i=1}^{N} (AT_{i}^{ref} - \overline{AT}^{ref})(AT_{i}^{c} - \overline{AT}^{c})}{\sqrt{\sum_{i=1}^{n} (AT_{i}^{ref} - \overline{AT}^{ref})^{2} \sum_{i=1}^{N} (AT_{i}^{c} - \overline{AT}^{c})^{2}}} \quad \text{and} \quad RE = \sqrt{\frac{\sum_{i=1}^{N} (AT_{i}^{ref} - \overline{AT}_{i}^{c})^{2}}{\sum_{i=1}^{N} (AT_{i}^{ref})^{2}}}$$

where $AT^{ref}$ and $AT^{c}$ denote respectively the reference AT and the computed ones. The number $N$ may be epicardial or endocardial nodes. The numbers $\overline{AT}^{ref}$ and $\overline{AT}^{c}$ are the mean values of $AT^{ref}$ and $AT^{c}$ respectively.

Correlation coefficients (CC) and Relative errors (RE) between the reconstructed AT and reference ones were computed. The AT obtained with both the VM and the MFS were compared to the reference, simulated, ones. In addition we could also compare the AT obtained with the VM to reference ones on the endocardium, including septum.

The robustness of the method was analyzed by adding different level of noise to the data. We added a gaussian noise to the data with different signal to noise ratio (SNR).

# 3 Results

In this section, we present the results obtained with the VM and compare them to classical MFS ones using these simulated datasets.

## 3.1 Data without noise

At first, we tested our VM described in section 2 and compared them with MFS on the datasets EDGAR. We report in Table 1 the correlation coefficient (CC) and the relative error (RE) of the activation times (AT) for both methods on all datasets. In addition, a boxplot of these results are depicted in figure 3. The figure 2 shows an example of the behavior of the reconstructed volume current source.

We see in Table 1 that the VM gives generally a better result in term of correlation coefficient (CC) and relative error (RE) than the MFS in these cases. Not that both the MFS and the VM are inaccurate on the SEPTCENTER dataset. For the data RVPOST, the MFS gives a very bad CC and RE. The best CC and RE are obtained for the data LVLAT for the VM (on the epicardium and the endocardium) while for the MFS it was for the data LVLATEPI. In figure 3, we observe that the VM is better than for the MFS (median CC : 0.75 for the VM and 0.4 for the MFS and median RE : 0.31 for the VM and 0.35 for the MFS) on the epicardium. In addition, the VM could reconstruct
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Example of the behavior of the reconstructed volume current source over time.</figcaption></figure>

<table><caption>**Table 1.** CC and RE on the AT for the MFS and the VM in all dataset (best result in bold font).</caption><thead><tr><th rowspan="3"></th><th colspan="4">VM</th><th colspan="4">MFS</th></tr><tr><th colspan="2">Epi</th><th colspan="2">Endo</th><th colspan="2">Epi</th><th colspan="2">Endo</th></tr><tr><th>CC</th><th>RE</th><th>CC</th><th>RE</th><th>CC</th><th>RE</th><th>CC</th><th>RE</th></tr></thead><tbody><tr><td>SEPTCENTER</td><td>0.26</td><td>0.32</td><td>0.36</td><td>0.40</td><td>0.17</td><td>0.40</td><td>-</td><td>-</td></tr><tr><td>LVLAT</td><td></td><td><strong>0.83</strong></td><td><strong>0.28</strong></td><td><strong>0.79</strong></td><td><strong>0.26</strong></td><td>0.69</td><td>0.31</td><td>-</td></tr><tr><td>LVAPEX</td><td>0.70</td><td>0.28</td><td>0.63</td><td>0.30</td><td>0.34</td><td>0.39</td><td>-</td><td>-</td></tr><tr><td>LVANT</td><td>0.74</td><td>0.29</td><td>0.66</td><td>0.36</td><td>0.47</td><td>0.33</td><td>-</td><td>-</td></tr><tr><td>RVPOST</td><td>0.75</td><td>0.37</td><td>0.37</td><td>0.54</td><td>-0.05</td><td>0.70</td><td>-</td><td>-</td></tr><tr><td>RVANT</td><td>0.67</td><td>0.37</td><td>0.70</td><td>0.38</td><td>0.10</td><td>0.38</td><td>-</td><td>-</td></tr><tr><td>LVLATEPI</td><td>0.78</td><td>0.35</td><td>0.73</td><td>0.28</td><td><strong>0.75</strong></td><td><strong>0.32</strong></td><td>-</td><td>-</td></tr><tr><td>LVLATENDO</td><td>0.79</td><td>0.31</td><td>0.66</td><td>0.29</td><td>0.70</td><td>0.30</td><td>-</td><td>-</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Boxplot of CC and RE on the AT for the VM and the MFS</figcaption></figure>

the activation times (AT) on the endocardium with acceptable CC and RE (median CC : 0.65 and median RE : 0.33). For the VM, the reconstruction on the epicardium is
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** AT maps for dataset LVLAT. A : reference. B : VM method. C : MFS method. Epicardium (left) and endocardium (right).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Scatter plot of AT reference versus reconstructed ones for LVLAT dataset. Left : Comparison between VM and MFS on the epicardium. Right : VM on the endocardium.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** AT maps for dataset SEPT. A : reference. B : VM method. C : MFS method. Epicardium (left) and endocardium (right).</figcaption></figure>

slightly better than on the endocardium (median CC : 0.75 vs 0.65 and median RE : 0.31 vs 0.33). For the datasets LVLAT and SEPTCENTER, we represent the AT maps (figure 4 and 6) and the scatter plot (figure 5 and 7). We can see that AT maps for the VM was smoother than for the MFS ones and we observe more clustering of points along horizontal lines for the MFS method. We remark also that late activation times
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Scatter plot of AT reference versus reconstructed ones for SEPT dataset. Left: Comparison between VM and MFS on the epicardium. Right: VM on the endocardium.</figcaption></figure>

were not well reconstructed for the dataset LVLAT (figure 5). For dataset SEPT, the scatter plot (figure 7) shows more dispersion that indicates the bad correlation between AT reference and reconstructed ones for both methods.

## 3.2 Data with noise

Next, we analyzed the sensitivity of both models with respect to noise by adding Gaussian noise with SNR from 5 to 20 dB to the original BSPM data. Figure 8 shows the mean of the correlation coefficient (CC) and the relative error (RE) obtained using different SNR on the epicardium with both methods. Figure 9 presents the results obtained on the endocardium for the VM. In figure 8, we observe a degradation of the correla-

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Mean of CC and RE with respect to the SNR for VM and MFS on the epicardium</figcaption></figure>

tion coefficient (CC) and the relative error (RE) for both methods on the epicardium when the SNR decreases. The CC (respectively the RE) vary from 0.27 to 0.60 (respectively 0.31 to 0.37) for the VM and from 0.05 to 0.24 (respectively 0.44 to 0.58) for the MFS. However, the VM is more robust than for MFS. We think that the classical torso formulation of the ECGI problem is more ill-posed than the formulation proposed.
<figure><img src="image_9.png" /><figcaption>Fig. 9. Mean of CC and RE with respect to the SNR for VM on the endocardium</figcaption></figure>

On the endocardium, we observe in figure 9 a similar behaviour for the VM. The CC (respectively RE) vary from 0.31 to 0.59 (respectively 0.31 to 0.37).

## 4 Discussion and conclusions

In this paper, we presented a volume method to solve inverse problem in electrocardiography. The proposed volume method reconstructs a whole-heart volume current source, instead of the epicardial potential field, as in the MFS. The VM is quite different from the MFS since it directly searches a volume current source instead of a potential, and also because it accounts for the electrical conductivity in the heart, and its possible variation in a scar, but also in the intracavitary blood, while the MFS only accounts for the torso conductivity.

For the datasets presented in section 3, the VM performed better than the MFS. In addition, the VM method could reconstruct acceptable endocardial activation (figure 3). Both methods degrade when the SNR decreases, but the VM always remains better than the MFS (figure 8).

All tests were completed with in silico data. It will be interesting to apply the VM to experimental or clinical data. It is also important to apply the VM method on more datasets with scars, as tested in a preliminary work [2]. At last, it would be interesting to study more regularization parameter chose for the VM method.

## Acknowledgments

This study received financial support from the French Government as part of the “Investments for the Future” program managed by the National Research Agency (ANR), Grant reference ANR-10-IAHU-04, and also from European Research Council (ERC), Grant ECSTATIC (715093).

## References

1. Ben Belgacem, F.: Why is the cauchy problem severely ill-posed? Inverse Problems **23**, 823 (03 2007). https://doi.org/10.1088/0266-5611/23/2/020
2. Diallo, M.M., Potse, M., Dubois, R., Coudière, Y.: Solving the ecgi problem with known locations of scar regions. In: Computing in Cardiology 2020, Rimini, Italy. https://doi.org/10.22489/CinC.2020.237

3. Duchateau, J., Potse, M., Dubois, R.: Spatially coherent activation maps for electrocardiographic imaging. IEEE Transactions on Biomedical Engineering 64(5), 1149–1156 (2017). https://doi.org/10.1109/TBME.2016.2593003

4. Hecht, F.: New development in freefem++. Journal of Numerical Mathematics 20(3-4), 251 – 266 (2012). https://doi.org/10.1515/jnum-2012-0013

5. Nielsen, B.F., Cai, X., Lysaker, M.: On the possibility for computing the transmembrane potential in the heart with a one shot method: An inverse problem. Mathematical Biosciences 210(2), 523 – 553 (2007). https://doi.org/10.1016/j.mbs.2007.06.003

6. Schulze, W.H., Potyagaylo, D., Schimpf, R., Papavassiliu, Tulumen, E., Rudic, B., Liebe, V., Doesch, C., Trachtler, J., Borggrefe, M., Dossel, O.: A simulation dataset for ecg imaging of paced beats with models for transmural, endo-and epicardial and pericardial source imaging. In: ECG Imaging 2015, Germany. https://doi.org/10.13140/RG.2.1.1946.8568

7. Wang, D., Kirby, R.M., MacLeod, R.S., Johnson, C.R.: Inverse electrocardiographic source localization of ischemia: An optimization framework and finite element solution. Journal of Computational Physics 250, 403 – 424 (2013). https://doi.org/10.1016/j.jcp.2013.05.027

8. Wang, Y., Rudy, Y.: Application of the method of fundamental solutions to potential-based inverse electrocardiography. Annals of Biomedical Engineering 34(8), 1272–1288 (2006). https://doi.org/10.1007/s10439-006-9131-7