# Comparison of Two Formulations for Computing Body Surface Potential Maps

Emma Lagracie¹, Lisl Weynans¹, Yves Coudière¹

¹ Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251, IHU Liryc, F-33400 Talence, France

## Abstract

In order to establish a link between the activation map and the body surface signals, we want to compute body potentials from activation maps and a frontlike approximation of the activation of the heart. To do so, two formulations that map the extracardiac potential from the transmembrane voltage naturally emerge from the bidomain model. Either the extracellular/extracardiac potential solves a Laplace equation with discontinuous conductivity coefficient and ionic current as a source (Source Formulation F1); or the quasi-stationary electrical balance between the intra- and extracellular fields (Balance Formulation F2). In this work, we compare F1 and F2, to determine which formulation is the most relevant to use.

We compute reference activation map $\psi$, transmembrane voltage $v$ and body surface map $u$ with a bidomain 2D code. We design two alternative shapes $\tilde{v}$ for a frontlike approximation of $v$. Afterwards, two extracellular/extracardiac potentials are computed from the activation time $\psi$ and $\tilde{v}$, using the two different formulations. Then the extracardiac potentials solutions of F1 and F2 are compared respectively to the solution of the bidomain $u$.

Results show that the Balance Formulation F2 is robust to the input data ($\tilde{v}$ and $\psi$). On the contrary, the Source Formulation F1 is very unstable and generates very large errors on the body surface map.

## 1. Introduction

In the framework of ECGi, we would like to reconstruct activation maps of the heart directly from torso potentials, available thanks to multi-electrode vests for example. To this end, we need to express a direct problem given the activation time $\psi(x)$. In particular, we want to deduce the extracellular and extracardiac potential $u$ directly from $\psi$. We choose to model the transmembrane voltage as an activation front with a predefined shape as in [1–3]. In this context, the approximate transmembrane voltage writes $\tilde{v}(x, t) = \tilde{V}(t - \psi(x))$. Even if its shape is uncertain, this representation of the transmembrane voltage captures the global behaviour of the depolarisation front.

Then, the extracellular and extracardiac potential is expressed as a solution of a partial differential equation with a source term that depends on $\tilde{V}$, $\psi$, and the formulation of the operator $[\tilde{v} \to u]$. This operator is derived from the bidomain equations that follow. In a heart ($\Omega_H$) embedded in a torso ($\Omega_T$) geometry, the bidomain model is given by

$$ 
\begin{aligned}
\operatorname{div}(\sigma_i \nabla u_i) &= \partial_t v + I_{\text{ion}}(v) && \text{in } \Omega_H, \\
\operatorname{div}(\sigma_e \nabla u) &= -\partial_t v - I_{\text{ion}}(v) && \text{in } \Omega_H, \\
\operatorname{div}(\sigma_T \nabla u) &= 0 && \text{in } \Omega_T, \\
\sigma_i \nabla (u + v) \cdot n &= 0 && \text{on } \partial \Omega_H, \\
\sigma_T \nabla u \cdot n &= 0 && \text{on } \partial \Omega_T,
\end{aligned}
\quad (1) $$

with $u_i \in H^1(\Omega_H)$, $u \in H^1(\Omega_H \cup \Omega_T)$ the intracellular and extracellular / extracardiac potentials respectively, $v = u_i - u \in \Omega_H$ the transmembrane voltage, and $\sigma_i$, $\sigma_e$, $\sigma_T$ the conductivities respectively in the cardiac intracellular medium, cardiac extracellular medium and torso medium. They are scaled to use nondimensionnal terms. The ionic current, $I_{\text{ion}}$ is a non-linear (often cubic) function of $v$ and other variables that depend on the ionic model chosen. Equivalently, the volumic equations in the heart can write

$$ 
\begin{aligned}
\operatorname{div}(\sigma_i \nabla (u + v)) &= \partial_t v + I_{\text{ion}}(v) && \text{in } \Omega_H, \\
\operatorname{div}((\sigma_i + \sigma_e) \nabla u) &= -\operatorname{div}(\sigma_i \nabla v) && \text{in } \Omega_H.
\end{aligned}
\quad (2) $$

This version is more widely used in the literature. In these two equivalent formulations of the bidomain we see two equations modelling the dependency between $u$ and $v$:

$$ 
\begin{aligned}
-\operatorname{div}(\sigma_e \nabla u) &= \partial_t v + f(v) && \text{in } \Omega_H, \\
-\operatorname{div}(\sigma_T \nabla u) &= 0 && \text{in } \Omega_T;
\end{aligned}
\quad (F1) $$

and

$$ 
\begin{aligned}
-\operatorname{div}(\sigma_e \nabla u) - \operatorname{div}(\sigma_i \nabla (u + v)) &= 0 && \text{in } \Omega_H, \\
-\operatorname{div}(\sigma_T \nabla u) &= 0 && \text{in } \Omega_T.
\end{aligned}
\quad (F2) $$

Thus, (F1) and (F2) can be seen as two independent problems which data is $v$ and that define two mappings of $[\tilde{v} \to u]$, taking $v = \tilde{v}$. Either $u$ solves a Laplace equation with discontinuous conductivity coefficients (heart
and torso) and a nonlinear source (Source Formulation (F1)), or it solves the quasi-stationary electrical balance between the intra and extracellular fields (Balance Formulation (F2)). Commonly, the potential $u$ is computed from the Balance Formulation (F2). Anyway, we may also use (F1), that has the advantage of not including any coupling on the heart boundary, and which differential operator is naturally shaped for the coupled heart/torso problem. Note that if $v$ solves the complete bidomain equations, the two formulations are equivalent. However, as soon as $v$ differs from the bidomain solution, (F1) and (F2) are no longer equivalent. In this study, we compare the solutions of the two formulations ($u_1$ for (F1) and $u_2$ for (F2)) evaluated in $\tilde{v}(x, t) = \bar{V}(t - \psi(x))$, with the bidomain solution (denoted by the index ref). We also investigate the robustness of those two formulations to small variations in the trans-membrane voltage and the activation time.

## 2. Method

### 2.1. Reference solution and activation map

Let ($u_{\text{ref}}$, $v_{\text{ref}}$) be the couple solution of the bidomain with the Mitchell-Shaeffer [4] model for the ionic potential, given by

$$ 
\begin{aligned}
\partial_t v &= \frac{1}{\tau_{in}}hv^2(1-v) - \frac{v}{\tau_{out}}, \\
\partial_t h &= \begin{cases} \frac{1-h}{\tau_{open}} & \text{is } v < v_{gate}, \\ \frac{-v}{\tau_{close}} & \text{otherwise.} \end{cases} \qquad (3)
\end{aligned}
$$

Concerning the numerical schemes, the 2D Discrete Duality Finite Volume (DDFV) [5] scheme was used for the space discretization and an order 2 Semi-implicit Backward Discretization scheme (SBDF2) for the temporal discretization. Thus, in the bidomain code, the term $\partial_t v + f(v)$ is approached by

$$ \frac{1}{dt} (1.5v^n - 2v^{n-1} + 0.5v^{n-2}) + 2f(v^{n-1}, h^{n-1}) - f(v^{n-2}, h^{n-2}), \quad (4) $$

for a time step $dt > 0$. Two 2D meshes were used, an artificial mesh of a circular heart in a circular torso, and a 2D slice of a segmentation of a patient's ventricle and torso.

The activation map of the heart was computed as a piecewise affine function given by its values $\psi_i$ at each node $x_i$ of the mesh (example Figure 3). The time $\psi_i(x_i)$ is defined as the first time at which $t \to v(x_i)$ crosses a threshold. This threshold is set at 0.5, as in the Mitchell-Shaeffer model $v$ takes values between 0 and 0.95. To find $\psi_i$ at each $x_i$, the functions $t \to v(x_i)$ are interpolated by cubic splines.

### 2.2. Numerical resolution of the two approached formulations

The equations (F1) and (F2) being static, we used again a 2D DDFV spatial scheme for the resolution and the discretization. Two predefined shapes $\bar{V}$ were designed. The first one is a Heaviside front, smoothed such that the right-hand sides (RHS) of the two formulations $\partial_t \tilde{v} + f(\tilde{v})$ and $-\text{div}(\sigma_i \nabla \tilde{v})$ are $L^2$ regular. The second one is the front obtained by solving the 0D Mitchell-Shaeffer problem, with a $C^\infty$ initial stimulation. For the formulation (F2), the RHS is computed by applying a DDFV scheme to the chosen form of $\tilde{v}$. For the source formulation (F1), the temporal derivative term is either explicitly given, $\partial_t \tilde{v} = \bar{V}'$, either derived numerically as in (4). In the Source Formulation (F1), the ionic term is a cubic function of $v$.

#### 2.2.1. Choice and computation of the front $\bar{V}$

In the first case, $\bar{V}$ is a smoothed Heaviside function depending on a front duration parameter $\epsilon$, specifically, $\tilde{v} = v_{\epsilon}(x, t) = \tilde{H}_{\epsilon}(t - \psi(x))$. To fulfill the constraint of a $H^1$ regularity for the solution of (F1), we imposed $\tilde{H}_{\epsilon}(x) = 0$ for $x < -\epsilon$, $\tilde{H}_{\epsilon}(x) = 0.94$ for $x > \epsilon$ and $\tilde{H}_{\epsilon}(x) = g$ in between ($|x| \le \epsilon$), with $g(x) = 0.94 (\frac{1}{4\epsilon^3}x^3 + \frac{3}{4\epsilon}x + \frac{1}{2})$. Both $v_{\epsilon}$ and its temporal derivative are then analytically known.

For the 0D Mitchell-Shaeffer front (MS0D), the Mitchell-Shaeffer model was solved in 0D, using classical python integration tools, on a time range of 330 ms (AP duration). Then, a cubic spline interpolation from python spycipy library was used to compute a model for $\tilde{v}$ callable at all times. The time scale was translated to center $\tilde{v}$ on the activation time.

#### 2.2.2. Computation of the ionic term

In the Mitchell-Shaeffer model, for an ideal propagation, $h = 1$ on the activation front. Then the Mitchell-Shaeffer function $f_{\text{MS}}$ becomes $f_{\text{MS}}(v) = \frac{v^2(1-v)}{\tau_{in}} - \frac{v}{\tau_{out}}$, with $\tau_{in} = 0.3$ and $\tau_{out} = 6$ as in the original paper [4]. Denoting $f_{\text{ref}}(v)$ the bidomain ionic term, we compared $f_{\text{MS}}$ to a cubic polynomial interpolation of $v_{\text{ref}} \to f_{\text{ref}}(v_{\text{ref}})$ denoted $f_{\text{int}}$. Several random points $x_k$ of the heart mesh were selected, and a least-square minimization was realised to find $v(x_k) \to f_{\text{int}}(v(x_k))$ under the constraint that 0 and 0.94 are zeros of the function. This means that two coefficients were left to find: the dominant coefficient and the third root. For all points $x_k$, the interpolated coefficients found were close to each other. Thus we took their mean value for designing the function $f_{\text{int}}(v)$. As illustrated on one mesh point in Figure 1 the $f_{\text{int}}$ function fits better the reference curve $f_{\text{ref}}(v_{\text{ref}})$ than $f_{\text{MS}}$ and was thus chosen for the following.
<figure><img src="image_1.png" /><figcaption>Figure 1. Comparison between $f_{MS}$ in green and $f_{int}$ in red, as functions of $v$. The stars are the reference values $f_{ref}(v)$ obtained by the bidomain propagation.</figcaption></figure>

# 3. Results

## 3.1. Verification of the resolution of the two formulations

To verify the resolution, we substituted $v_{ref}$ in both formulations. For the formulation (F1), the RHS is computed with (4). $u_1$ and $u_2$ are expected to match exactly the reference $u_{ref}$. In practice, the differences are $\approx 10^{-13}$ due to numerical errors. However, the result $u_1$ is very sensitive to small errors in the RHS. For instance, if we replace formula (4) by another discretization scheme of same order, the errors $\frac{\|u_1 - u_{ref}\|_{L^2(\Omega_T)}}{\|u_{ref}\|_{L^2(\Omega_T)}}$ jump to an order $10^{-2}$ or even $10^{-1}$, as presented in Tables 1 and 2 for the circular mesh.

<table><caption>Table 1. Relative $L^2(\Omega_T)$ errors for the solution of (F1) when changing $f$ in (4).</caption><thead><tr><th>Ionic term</th><th><div style="text-align: center;"><span style="display: block; border-bottom: 1px solid black; padding-bottom: 1px;">||u<sub>1</sub> - u<sub>ref</sub>||<sub>L<sup>2</sup>(&Omega;<sub>T</sub>)</sub></span><span style="display: block; padding-top: 1px;">||u<sub>ref</sub>||<sub>L<sup>2</sup>(&Omega;<sub>T</sub>)</sub></span></div></th></tr></thead><tbody><tr><td>f<sub>ref</sub> (v<sub>ref</sub>) (reference)</td><td>1.5E-13</td></tr><tr><td>f<sub>int</sub> (v<sub>ref</sub>(t<sup>n</sup>))</td><td>8.1E-01</td></tr><tr><td>f<sub>MS</sub> (h(t<sup>n</sup>), v<sub>ref</sub>(t<sup>n</sup>))</td><td>4.4E-02</td></tr></tbody></table>

These results anticipated mediocre performances for (F1).

<table><caption>Table 2. Relative $L^2(\Omega_T)$ errors for the solution of (F1) when changing $\partial_t v$ in (4).</caption><thead><tr><th>Discrete derivative scheme</th><th><div style="text-align: center;"><span style="display: block; border-bottom: 1px solid black; padding-bottom: 1px;">||u<sub>1</sub> - u<sub>ref</sub>||<sub>L<sup>2</sup>(Ω<sub>T</sub>)</sub></span><span style="display: block; padding-top: 1px;">||u<sub>ref</sub>||<sub>L<sup>2</sup>(Ω<sub>T</sub>)</sub></span></div></th></tr></thead><tbody><tr><td>SBDF2 (reference)</td><td>1.5E-13</td></tr><tr><td>Euler centered</td><td>2.5E-02</td></tr><tr><td>Explicit Euler</td><td>4.3E-01</td></tr></tbody></table>

in case of a pre-shaped front, the source term being very

<figure><img src="image_2.png" /><figcaption>Figure 2. Bidomain solution $u_{ref}$ on the left, and $u_1$ solution of (F1) on the right for $\partial_t v = \partial_t v_{ref}$ and $f(v) = f_{int}(v_{ref}(t^n))$ in (4) (line 2 of Table 1).</figcaption></figure>

sensitive to small variations, and relying on a balance between the transmembrane voltage time derivative and the ionic term. Figure 2 illustrates this observation.

## 3.2. Analysis of the effects of the formulation

<figure><img src="image_3.png" /><figcaption>Figure 3. Activation times computed from the bidomain propagation.</figcaption></figure>

In the following, $v$ is substituted by $\tilde{v} = v_{\epsilon}$ or the MS0D front. The relative differences between $u_1$ and $u_{ref}$ and $u_2$ and $u_{ref}$ are presented in Figure 4. Accordingly to the previous results, (F1) generates huge errors, thus cannot be used in computations. On the contrary, (F2) seems robust and generates less than 10% of difference for all the tested $\tilde{v}$. For parameter $\epsilon$ between 0.5 and 5 for $v_{\epsilon}$, the errors are lower than the ones obtained with the MS0D model. We can also notice that there exist an optimal $\epsilon_0$ parameter that minimizes the difference $\frac{\|u_2 - u_{ref}\|_{L^1(\partial\Omega_T \times (0, T))}}{\|u_{ref}\|_{L^1(\partial\Omega_T \times (0, T))}}$. Anyway, if we look at the difference in space only, $\epsilon_0$ slightly depends on time.

## 3.3. First sensitivity analysis

Focusing on (F2), a white noise perturbation $w$ was added to the RHS of the equation. Four noisy trans-
<figure><img src="image_4.png" /><figcaption>Figure 4. Relative differences between $u_1$ (red), $u_2$ (blue) and $u_{\text{ref}}$ respectively, and $\tilde{v}$ and $v_{\text{ref}}$ (black). $\tilde{v}$ may be the MS0D model (horizontal dashed lines), or $v_{\epsilon}$ with varying duration parameter $\epsilon$ (points).</figcaption></figure>

membrane voltages were designed: $\tilde{v} = v_{\text{ref}} + w(x)$, $\tilde{v} = v_{\epsilon} + w(x)$, $\tilde{v} = H_{\epsilon}(t - \psi(x) + w(x))$ and $\tilde{v} = H_{\epsilon}(t - \psi(x) + w_0)$, where $w_0$ is a scalar. The last case represent a displacement of the whole activation front line, as it would be if the initial activation site were uncertain. For the % of noise on the activation time, a reference amplitude was needed. Denoting $p$ the position of the stimulation site, we allow $p + \Delta p$ to be located in all septum. Then $\Delta t = \frac{\Delta p}{c}$ with $c$ the mean speed of the activation front in the septum. Surprisingly, the chart Figure 5 shows

<figure><img src="image_5.png" /><figcaption>Figure 5. Relative difference between $u_2$ and $u_{\text{ref}}$ on the torso at time $t = 35$ ms for 5000 realisations of white noise $w$. $\epsilon = 2.5$.</figcaption></figure>

no significant difference between the solutions of (F2) with $\tilde{v} = v_{\text{ref}} + w$ and $\tilde{v} = v_{\epsilon} + w$, confirming again the validity of using (F2) in future computations. Moreover, the

formulation does not appear to be very sensitive to small noise in the stimulation site. In this case, the errors are significantly inferior to raw gaussian noise on the activation time.

## 4. Discussion and Conclusion

Our results point out clearly that the Source Formula- tion (F1) cannot be used for propagating $u$ from $v$. Even with very accurate entries, the body surface potential map is no longer preserved. This emphasizes the underlying equilibrium existing between the ionic term and $v$. On the contrary, the Balance Formulation allows to correctly recover body surface potentials even with a simple preshaped transmembrane potential or with noise.

It would be interesting to push further this study in 3D. Moreover, a real sensitivity analysis of $u_2$ as function of $\psi$ would be very useful. In future work, we could also study the sensitivity of formulation (F2) to heterogeneities in the heart conductivity.

## Acknowledgments

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

## References

[1] Franzone PC, Guerri L, Tentoni S. Mathematical modeling of the excitation process in myocardial tissue: Influence of fiber rotation on wavefront propagation and potential field. Mathematical biosciences 1990;101(2):155-235.

[2] Van Dam PM, Oostendorp TF, Linnenbank AC, Van Oosterom A. Non-invasive imaging of cardiac activation and recovery. Annals of biomedical engineering 2009;37:1739-1756.

[3] Ravon G, Coudière Y, Potse M, Dubois R. Impact of the endocardium in a parameter optimization to solve the inverse problem of electrocardiography. Frontiers in Physiology 2019;9:1946.

[4] Mitchell CC, Schaeffer DG. A two-current model for the dynamics of cardiac membrane. Bulletin of mathematical biology 2003;65(5):767-793.

[5] Delcourte S, Domelevo K, Omnes P. Discrete-duality finite volume method for second order elliptic problems. HAL 2005;2005.

Address for correspondence:

Emma Lagracie
200 avenue de la vieille tour, 33400 Talence, France
emma.lagracie@inria.fr