# An Epicardial Model with a Physiological Total Variation Regularization for Solving the ECGi Problem

Emma Lagracie, Yves Coudière, Lisl Weynans

<figure><img src="image_1.png" /></figure>

## To cite this version:

Emma Lagracie, Yves Coudière, Lisl Weynans. An Epicardial Model with a Physiological Total Variation Regularization for Solving the ECGi Problem. 2026. hal-05521199

HAL Id: hal-05521199

https://hal.science/hal-05521199v1

Preprint submitted on 20 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An Epicardial Model with a Physiological Total Variation Regularization for Solving the ECGi Problem

Emma Lagracie¹ and Yves Coudière²,³ and Lisl Weynans²,³

¹ Otto von Guericke Universität, Magdeburg, Germany

² Univ. Bordeaux, IMB, UMR 5251, IHU Liryc, F-33400 Talence, France

³ CARMEN team, Inria Bordeaux, F-33400 Talence, France

February 20, 2026

## Abstract

Electrocardiographic imaging (ECGi) aims to reconstruct cardiac electrical activity from noninvasive measurements of torso surface potentials. This inverse problem is commonly formulated as a PDE-constrained optimization problem relying on a source model. One common way to overcome its severe ill-posedness is to add a regularization term to the cost function to be minimized in the optimization problem. Another way consists in adding constraints in the optimization problem, for instance based on physical informations or modeling. Following this idea, in this work, we introduce a novel source model, named the epicardial model (EPI), which provides a heuristic three-dimensional extension of a previously proposed depth-averaged two-dimensional model. This new model is formulated in the torso domain only and couples a bidomain-like surface equation on the heart with the standard Laplace equation in the torso. Thus, it can naturally be substituted into standard torso-based clinical workflows while offering more possibilities for incorporating physiologically relevant prior information on the modeled epicardial surface.

The epicardial model is first validated numerically against reference bidomain simulations. While amplitude discrepancies are observed, the model is able to reproduce the spatial distribution of the electrical potentials. The epicardial model is then evaluated within an ECGi framework. Two EPI-constrained optimization formulations are considered: one using classical Tikhonov regularization and another using a physiological total variation (TV) regularization applied to the transmembrane voltage.

Reconstructions are compared with those obtained from the classical Laplace-based Cauchy formulation. Without addition of information, the proposed EPI based framework yields reconstruction of comparable accuracy to standard approaches despite the presence of modeling approximations. The incorporation of a prior through the TV regularization, allowed by the model, improves signals reconstruction and reduces artificial lines of block in activation maps. These results highlight the potential of the epicardial model as a physiologically enriched alternative to classical ECGi formulations.

**Keywords** ECGi, Inverse problem, Source model, Total variation regularization, Cardiac electrophysiology, Epicardial model

## 1 Introduction

Electrocardiographic imaging, also called ECGi, aims at non-invasively recovering cardiac bioelectrical signals from measured body surface potentials. In particular, the main goal is often to reconstruct the electrical activation sequence of the heart, which may provide insight on the cause of pathological cardiac rhythm. However, the ECGi problem is mathematically ill-posed in the sense of Hadamard, meaning that at least one of the following conditions fail: existence of a solution, uniqueness of a solution, continuous dependency to the data. In practice, we observe high instability in the solutions with respect to the data, which makes this inverse problem particularly difficult to solve.

Usually, ECGi is formulated as a partial differential equation constrained optimization problem [2, 22, 27]. Subject to the "source model" establishing a link between bioelectric sources within the heart (eg. transmembrane voltage or extracellular potentials) and observable torso surface potentials, one seeks to find the cardiac signals producing body surface potentials as close as possible
to the measured data [19, 4]. A regularization term is added, ensuring the well-posedness of the minimization problem. Most often, static reconstructions corresponding to different time instants are calculated, then post-processed to compute an activation map representing the activation sequence of the heart [10, 30].

Source models as in references [5, 12, 26] can be derived from the bidomain equations [32, 13], which constitute the reference equations in cardiac electrophysiology macroscopic modeling. They can involve a volume or surface heart domain associated with different corresponding bioelectric signals (transmembrane voltage and/or extracellular potential), coupled with an extracardiac potential defined in the torso domain. Usual source models are the electrostatic equilibrium of the volume bidomain model [22, 9, 34, 16, 38], the surface equivalent dipole layer (EDL) model [12, 33, 36, 41], and the Laplace equation in the torso [2, 4], which uses only a surface representation of the heart.

The choice of the source model used in the optimization problem has a real impact on the nature of the inverse problem, and can change the reconstructions. For instance, the use of the electrostatic equilibrium equations in ECGi allows to recover both extracellular potential and transmembrane voltage within the whole myocardium volume, while incorporating a priori information on the heart electrical conductivities. However, this formulation also implies a lack of uniqueness of the solutions, in addition to the lack of existence common to all ECGi formulations [18, 17]. Surface formulations, like the Equivalent Dipole Layer formulation or the Laplace equation formulation, guarantee the uniqueness of a solution up to a constant, but do not offer to recover both the extracellular potential $u$ and the transmembrane voltage $v$ on the heart surface. They also leave less room for incorporating a priori information on tissue conductivities, or complex regularization patterns on $v$ and $u$ for instance. In practice, the most widely used formulation of ECGi, already applied in clinics consists in solving the static Cauchy problem based on the Laplace equation in the torso to recover the epicardial potential distribution on the heart surface [3, 7, 28, 39, 23, 6].

In this paper, we propose an alternative surface source model to be used as constraint in the inverse problem of ECGi. The proposed model heuristically extends to three dimensions the two-dimensional formulation developed in a previous work [21]. It combines a volume description of the potential field in the torso with a surface description of the transmembrane voltage and extracellular potential on the epicardium. Heart-torso transmission effects are incorporated through a Robin boundary condition. The so called epicardial source model can be interpreted as a static version of the bidomain model averaged over a thin myocardial layer.

Compared with volumetric approaches [9, 34], the proposed formulation reduces the number of unknowns to be reconstructed and alleviates the non-uniqueness issues associated with transmembrane voltage solutions. Within a heart surface-based framework, it enables the manipulation of bidomain-like equations and their associated modeling capabilities while relying solely on a torso mesh. This modeling framework therefore extends classical Laplace-based models by allowing the incorporation of physiologically relevant features, while preserving the practical advantage of using the same types of data and mesh. As a result, it can readily substitute the classical Laplace-based formulation in clinical oriented technologies.

Our epicardial model is primarily numerically validated against electrical signals simulated using the bidomain model, which serves as the reference. Then an epicardial model-based ECGi framework is formulated as a PDE-constrained optimization problem supplemented with a regularization term acting on the transmembrane voltage source. In this work, we investigate both a classical first-order Tikhonov regularization and a total variation (TV) regularization, which is made possible by the explicit representation of the transmembrane voltage in the epicardial model. For the Tikhonov-regularized problem, the optimality system is derived within a Lagrangian framework, leading to the solution of a linear system in practice. The TV-regularized formulation is solved using an iterative algorithm involving the resolution of Lagrangian subsystems [35, 14].

The introduction of TV regularization on the transmembrane voltage is physiologically motivated and illustrates one of the advantages of the epicardial model, namely the ability to incorporate meaningful prior information at the source level. Over the past decade, the ECGi inverse problem community has shown increasing interest in TV regularization [14, 37, 40, 31, 15], which consists in minimizing the $L^1$ norm of the gradient of the source term. Unlike $H^1$-type regularization, TV regularization promotes piecewise-constant solutions with sharp localized variations. Previous studies have mainly applied this strategy to classical potential-based formulations, with the goal of reconstructing extracardiac potentials exhibiting more physiologically realistic features. While the extracellular potential may exhibit steep variations at the activation front combined with smoother spatial variations elsewhere, the transmembrane voltage typically presents a plateau-like spatial profile during the activation phase of the heart, making it especially well suited to TV regularization.

The remainder of the article is organized as follows. Section 2 presents the reference bidomain for-
ward model and classical ECGi source models, and then introduces the proposed three-dimensional epicardial model together with its natural finite element discretization and the corresponding inverse ECGi formulation. In particular, three alternative surface-based ECGi formulations are considered: the classical Laplace equation-based model with Tikhonov regularization, the proposed epicardial model with Tikhonov regularization, and the epicardial model combined with total variation regularization.

Section 3 describes the methodology used to evaluate both the forward and inverse formulations of the epicardial model. It details the forward and inverse computational workflows, the parameters used to solve the bidomain model employed to generate reference signals, the procedure for constructing data for the inverse problem, and the evaluation metrics used to assess model performance. This section also specifies the chosen parameters of the epicardial model, including electrical conductivities and tunable coefficients, as well as the different scenarios considered for evaluating the model in both forward and inverse settings.

In Section 4, we numerically evaluate the epicardial forward model, and assess its use in an inverse setting. Finally, on Section 5 epicardial model-based ECGi reconstructions are assessed and compared with standard Laplace-Tikhonov ECGi solutions, for simulated data. The article concludes with a discussion in Section 6.

## 2 Forward and inverse models

We consider an open domain $\Omega \subset \mathbb{R}^3$, composed of three disjoint open subsets: the blood cavities $\Omega_B$, the cardiac muscle $\Omega_H$, and the torso $\Omega_T$, such that $\Omega_B \subset \Omega_H$ (meaning that $\bar{\Omega}_B \subset \Omega_H$), $\Omega_H \subset \Omega$, and $\bar{\Omega} = \bar{\Omega}_B \cup \bar{\Omega}_H \cup \bar{\Omega}_T$. The interface between the cardiac muscle and the torso is called the epicardium, denoted by $\Gamma_H = \bar{\Omega}_H \cap \bar{\Omega}_T$, while the interface between the blood and the cardiac muscle is called the endocardium. The torso surface, called body surface is the boundary $\Gamma_T = \partial\Omega$, see Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: Blood $\Omega_B$, cardiac $\Omega_H$, and torso $\Omega_T$ subdomains, epicardium $\Gamma_H$, and body $\Gamma_T$ surfaces, and normal vectors $n_H$ and $n_T$.</figcaption></figure>

### 2.1 Bidomain equations

The standard system of partial differential equations used to model the propagation of cardiac action potential, and the extracardiac electrical field is the bidomain equations:

$$ \operatorname{div}(\sigma_i \nabla u_i) = \chi (c_m \partial_t v + I_{\text{ion}}(v, w) + s) \quad \text{in } \Omega_H, \quad (1a) $$

$$ -\operatorname{div}(\sigma_e \nabla u) = \chi (c_m \partial_t v + I_{\text{ion}}(v, w) + s) \quad \text{in } \Omega_H, \quad (1b) $$

$$ -\operatorname{div}(\sigma_T \nabla u) = 0 \quad \text{in } \Omega_T \cup \Omega_B, \quad (1c) $$

where we look for the electrical fields $u_i(t) \in H^1(\Omega_H)$ and $u(t) \in H^1(\Omega)$. The function $v(t) = u(t) - u_i(t) \in H^1(\Omega_H)$ is called the transmembrane voltage. The function $w(t, x) \in \mathbb{R}^m$, defined for $x \in \Omega_H$ contains the state variables of the cardiac tissue, and the function $I_{\text{ion}}: \mathbb{R} \times \mathbb{R}^m \to \mathbb{R}$ is called the ionic current. $\chi$ is the surface to volume ratio of cell membranes and $c_m$ is the cell membrane surfacic capacitance. $\sigma_i$ is the tensorial intracellular conductivity, and $\sigma_e$ the tensorial extracellular conductivity. The source term $s(t, x)$ is an externally applied current, used only temporarily to trigger a local action potential. The state variables solve a system of ordinary differential equations called the ionic model. Here, we consider the simple ($m = 1$) Mitchell-Schaeffer model (MS) [24]:

$$ I_{\text{ion}}(v, w) = \frac{1}{\tau_{\text{in}}} w v^2 (1 - v) - \frac{v}{\tau_{\text{out}}}, \quad \partial_t w = \frac{h_{\infty}(v, w) - h}{\tau(v, w)}, \quad (2) $$
and $h_\infty(v, w) = 1$ if $v < v_g$ and 0 otherwise, $\tau_\infty(v, w) = \tau_{\text{open}}$ if $v < v_g$ and $\tau_{\text{close}}$ otherwise. The equations are completed by homogeneous boundary conditions

$$ \sigma_i \nabla u_i \cdot n_H = 0 \quad \text{on } \Gamma_H \cup \partial \Omega_B, \quad \sigma_T \nabla u \cdot n_T = 0 \quad \text{on } \Gamma_T, \qquad (3) $$

and the initial condition $v(0) = v_0$ and $w(0) = 0$, in $\Omega_H$. These equations have a solution (see e.g. [5]), and the trace of $u$ on $\Gamma_T$ yields the well known electrocardiogram. They are used to solve direct problems, i.e. to compute propagation of action potentials, as triggered by some local (in time and space) applied current $s$.

We compute approximate solutions of (1) using the standard Lagrange P1 finite element method on a simplicial mesh of $\Omega$, conformal with the subdomain $\Omega_H$ and in $\Omega$. It yields:

$$ M (c_m V' + I_{\text{ion}}(V, W) + S) + K_i U_i = 0 \qquad (4a) $$

$$ - \begin{pmatrix} M \\ 0 \end{pmatrix} (c_m V' + I_{\text{ion}}(V, W) + S) + K U = 0. \qquad (4b) $$

The matrices $M$, $K_i$, and $K$ are the usual mass and stiffness matrices, in $\Omega_H$, and $\Omega_B \cup \Omega_T$, with entries given by:

$$ [M]_{ij} = \int_{\Omega_H} \chi \phi_j \phi_i, \quad [K_i]_{ij} = \int_{\Omega_H} \sigma_i \nabla \phi_j \cdot \nabla \phi_i, \quad 1 \le i, j \le N_H, \\ [K]_{ij} = \int_{\Omega} \sigma \nabla \phi_j \cdot \nabla \phi_i, \quad 1 \le i, j \le N. $$

Here $N_H$ is the number of mesh points in $\Omega_H$, and $N$ is the number of mesh points in $\Omega$. The conductivity coefficient $\sigma$ is defined piecewise on the subdomains, $\sigma = \sigma_e$ in $\Omega_H$, $\sigma = \sigma_B$ in $\Omega_B$, and $\sigma = \sigma_T$ in $\Omega_T$. The functions $\phi_i$ are the basis functions of the P1 Lagrange space, the vector of unknowns are $U_i \in \mathbb{R}^{N_H}$ and $U \in \mathbb{R}^N$ (with $V = U_i - U_{|\Omega_H}$). We finally use the S-BDF2 time-stepping method proposed in reference [11] in order to compute time sequences of solutions $(V^n, U^n)$.

## 2.2 Source models for the ECGi problem

### 2.2.1 Bidomain electrostatic source model

Remark that $u \in H^1(\Omega)$ and $u_i \in H^1(\Omega_H)$ which solve equations (1) and boundary conditions (3) are such that $\sigma_i \nabla(u + v) \cdot n_H = 0$, and $\sigma_e \nabla u \cdot n_H = \sigma_T \nabla u \cdot n_H$ on $\Gamma_H$, and $u$ is continuous through $\Gamma_H$. We can compute the difference between equations (1b) and (1a), and their associated boundary conditions on $u$ and $u_i = u + v$. Given a transmembrane voltage $v \in H^1(\Omega_H)$, we observe that $u \in H^1(\Omega)$ solves the following quasi-static electrical balance equation:

$$ - \operatorname{div} ((\sigma_e + \sigma_i) \nabla u) - \operatorname{div} (\sigma_i \nabla v) = 0 \qquad \text{in } \Omega_H, \qquad (5a) $$

$$ - \operatorname{div} (\sigma \nabla u) = 0 \qquad \text{in } \Omega_T \cup \Omega_B, \qquad (5b) $$

with the transmission conditions

$$ (\sigma_i + \sigma_e) \nabla u \cdot n_H + \sigma_i \nabla v \cdot n_H = \sigma_T \nabla u \cdot n_H \quad \text{on } \Gamma_H \cup \partial \Omega_B. \qquad (6) $$

Given $v \in H^1(\Omega_H)$, by the Lax-Milgram theorem, equation (5) with the conditions (6) and $\sigma_T \nabla u \cdot n_T = 0$ has a unique solution, up to an additive constant.

It differs from the complete bidomain equations because here $v$ is an arbitrary function of $H^1(\Omega_H)$, meaning that $u$ and $v$ may not verify the two boundary conditions $\sigma_i \nabla(u + v) \cdot n_H = 0$, and $\sigma_e \nabla u \cdot n_H = \sigma_T \nabla u \cdot n_H$ separately on $\Gamma_H$, but only (6) instead.

Equations (5) with the conditions (6) and $\sigma_T \nabla u \cdot n_T = 0$ on $\Gamma_T$ constitute a three-dimensional source model for the ECGi problem, in which the source is the transmembrane voltage $v$, see reference [9] for an example application.

For this static equation, the P1 Lagrange discrete problem reads:

$$ \begin{pmatrix} K_i & 0 \\ 0 & 0 \end{pmatrix} U + K U = - \begin{pmatrix} K_i \\ 0 \end{pmatrix} V, \qquad (7) $$

where the matrices $K_i$ and $K$ are the same as for equations (4). In this case, the vector of unknowns is $U \in \mathbb{R}^N$, and the data is $V \in \mathbb{R}^{N_H}$. Usually a small mass matrix is added to this equation to impose unicity of the linear system to solve.
### 2.2.2 The standard torso only source model (STD)

In addition, we recall one of the most used [3, 7] source model, which consist simply of (1c), or equivalently (5b) with the no-flux boundary condition $\sigma_T \nabla u \cdot n_T = 0$ on $\Gamma_T$, and the Dirichlet condition $u = u_H \in L^2(\Gamma_H)$ on the epicardium.

$$ - \operatorname{div} (\sigma_T \nabla u) = 0 \qquad \text{in } \Omega_T, \qquad (8a) $$

$$ u = u_H \qquad \text{on } \Gamma_H, \qquad (8b) $$

$$ \sigma_T \partial_n u = 0 \qquad \text{on } \Gamma_T. \qquad (8c) $$

This system has obviously a unique solution $u \in H^1(\Omega_T)$. In the discrete setup, it yields the linear system

$$ \tilde{K}_T U = \begin{pmatrix} U_H \\ 0 \end{pmatrix}, \qquad (9) $$

where the matrix $\tilde{K}_T$ is the matrix $K_T$ defined as follows modified to account for the Dirichlet boundary condition $u = u_H$ on $\Gamma_H$.

$$ [K_T]_{ij} = \int_{\Omega_T} \sigma \nabla \phi_j \cdot \nabla \phi_i, \quad 1 \le i, j \le N_T. $$

$U_H \in \mathbb{R}^{N_{\Gamma_H}}$ is the vector of the degrees of freedom of $u_H$. Here $N_T$ is the number of mesh points in $\Omega_T$. The unknowns is the vector $U \in \mathbb{R}^{N_T}$. Note that $U$ is the extracardiac potential defined only on the torso subdomain $\Omega_T$, unlike in section 2.2.1 where it is defined on the whole $\Omega$.

### 2.2.3 A new epicardial source model (EPI)

It is usually observed on body surface electrical signals that variations due to activation of the cardiac muscle can be observed only after a breakthrough of the action potential through the epicardial surface, meaning that endocardial activation is below the observation range.

Hence, we propose a new source model, that can be interpreted as follows: the cardiac muscle is limited to a certain observable depth $h(x)$ ($x \in \Gamma_H$) from the epicardial surface, defined as $h(x) = p \times \text{wall thickness at } x$ for $0 < p < 1$; the fibers remain (almost) parallel to $\Gamma_H$ within this depth; we can average the transmembrane voltage $v$, and the electrical potential $u$ along the depth, but keep the electrical potential in $\Omega_T$, see Figure 2. Following the approach developed in [21], this

<figure><img src="image_4.png" /><figcaption>Figure 2: Schematic representation of the different elements appearing in the epicardial model.</figcaption></figure>

leads to the following equations and transmission condition for the unknown $\bar{u} \in H^1(\Gamma_H)$, called the effective electrical potential, $\bar{v} \in H^1(\Gamma_H)$, called the effective transmembrane voltage, the electrical potential $u \in H^1(\Omega_T)$ in the torso:

$$ \operatorname{div}_s (h(\bar{\sigma}_i + \bar{\sigma}_e) \nabla_s \bar{u}) - \sigma_T \nabla u \cdot n_H = - \operatorname{div}_s (h(\bar{\sigma}_i \nabla_s \bar{v})) \qquad \text{in } \Gamma_H, \qquad (10a) $$

$$ \operatorname{div} (\sigma_T \nabla u) = 0 \qquad \text{in } \Omega_T, \qquad (10b) $$

$$ \sigma^n \frac{u - \bar{u}}{\alpha h} = - \sigma_T \nabla u \cdot n_H \qquad \text{on } \Gamma_H, \qquad (10c) $$

where the normal vector $n_H$ is outward of $\Omega_T$. Equation (10a) is obtained by integration of equation (5a) along the depth $h(x)$ ($x \in \Gamma_H$), assuming no flux of current on the endocardial side, and using the transmission condition (6); equation (10b) is unchanged from (5b), and equation (10c) is a closure equation in which we approximate the flux of current at the epicardium as the finite difference $\sigma^n \frac{u - \bar{u}}{\alpha h}$. Of course, the no-flux boundary condition $\sigma_T \nabla u \cdot n_T = 0$ remains on $\Gamma_T$.
The coefficient $\sigma^n$ is defined as $\sigma^n = n_H \cdot (\sigma_i + \sigma_e)n_H$ and $\alpha$ satisfies $0 < \alpha < 1$. The operators $\text{div}_s$, and $\nabla_s$ denote the divergence and gradient operators on the surface $\Gamma_H$. The $2 \times 2$ matrices $\bar{\sigma}_i(x)$ and $\bar{\sigma}_e(x)$ are effective surface intracellular and extracellular conductivity tensors. Ideally, they should be such that we can write at any $x \in \Gamma$, by integration along the depth direction $z$, $\frac{1}{h} \int_0^h (\sigma_i + \sigma_e) \nabla u dz = (\bar{\sigma}_i + \bar{\sigma}_e) \nabla_s \bar{u}$, and $\frac{1}{h} \int_0^h \sigma_i \nabla v dz = \bar{\sigma}_i \nabla_s \bar{v}$. A natural choice is depth averages of the three-dimensional bidomain matrices, $\sigma_i$ and $\sigma_e$, although we have observed that it may not be the best choice [20].

Given $\bar{v} \in H^1(\Gamma_H)$, equations (10) define a unique solution $(\bar{u}, u) \in V := (H^1(\Gamma_H) \times H^1(\Omega_T)) / \mathbb{R}$ (i.e. $\bar{u}$ and $u$ are defined up to the same additive constant) by the Lax-Milgram theorem, based on the weak formulation:

$$ (\bar{\phi}, \phi) \in V, \quad \int_{\Gamma_H} h(\sigma_i + \bar{\sigma}_e) \nabla_s \bar{u} \nabla_s \bar{\phi} + \int_{\Omega_T} \sigma_T \nabla u \cdot \nabla \phi + \int_{\Gamma_H} \frac{\sigma^n}{\alpha h} (\bar{u} - u) (\bar{\phi} - \phi) \\ = - \int_{\Gamma_H} h \bar{\sigma}_i \nabla_s \bar{v} \nabla_s \bar{\phi}. \quad (11) $$

Equations (10) with the condition $\sigma_T \nabla u \cdot n_T = 0$ on $\Gamma_T$ constitute the epicardial source model for the ECGi problem, in which the source is the epicardial transmembrane voltage $\bar{v}$, see reference [21] for a two-dimensional version.

This model couples the epicardial surface to the torso while neglecting fluxes originating from the intracardiac cavities, and introduces effective surface conductivities as well as a variable thickness parameter. Fully substitutable for the classical approach based on the Laplace equation in the torso, this source model no longer allows reconstruction of transmural information, but it retains the ability to jointly reconstruct the extracellular potential $u$ and the transmembrane voltage $v$ on the epicardium.

For this static equation, the P1 Lagrange discrete problem reads:

$$ (\bar{K}_i + \bar{K}_e)\bar{U} + \frac{\sigma^n}{\alpha}\bar{M}(\bar{U} - \Pi_{\Gamma_H}U) = -\bar{K}_i\bar{V} \quad (12a) $$

$$ K_T U - \left( \begin{array}{c} \frac{\sigma^n}{\alpha} \\ 0 \end{array} \right) (\bar{U} - \Pi_{\Gamma_H} U) = 0. \quad (12b) $$

where the matrices $\bar{K}_{i,e}$, $K_T$, and $\bar{M}$ are again stiffness and mass matrices, but now $\bar{K}_{i,e}$ and $\bar{M}$ are defined on the discrete surface $\Gamma_H$, and $K_T$ is defined in the torso volume $\Omega_T$ only:

$$ [\bar{M}]_{ij} = \int_{\Gamma_H} \frac{1}{h} \phi_j \phi_i, \quad [\bar{K}_{\bullet}]_{ij} = \int_{\Gamma_H} h \bar{\sigma}_{\bullet} \nabla \phi_j \cdot \nabla \phi_i, \quad 1 \le i, j \le N_{\Gamma_H}, \bullet \in \{i, e\}, \\ [K_T]_{ij} = \int_{\Omega_T} \sigma \nabla \phi_j \cdot \nabla \phi_i, \quad 1 \le i, j \le N_T. $$

Here $N_{\Gamma_H}$ is the number of mesh points in $\Gamma_H$, and $N_T$ is the number of mesh points in $\Omega_T$. The projection matrix $\Pi_{\Gamma_H}$ is the discrete trace operator from $\Omega_T$ to $\Gamma_H$: $\Pi_{\Gamma_H,ij} = 0$, for $1 \le i \le N_{\Gamma_H}$, $1 \le j \le N_T$, except if $i$ is such that node $i$ belongs to $\Gamma_H$, in which case $\Pi_{\Gamma_H,ii} = 1$. The unknowns are the vectors $\bar{U} \in \mathbb{R}^{N_{\Gamma_H}}$, $\bar{U} \in \mathbb{R}^{N_T}$. For this problem, the vector $\bar{V} \in \mathbb{R}^{N_{\Gamma_H}}$ is a given data and not an unknown to compute.

To ensure the invertibility of the system, we use the usual technique of adding a scaled mass matrix on $\Gamma_H$ in the equation.

## 2.3 Formulations of Inverse problem for ECGi

In this section we look for electrical information in the cardiac muscle from the knowledge of the epicardial and body surfaces, and time series of body surface electrical potential data. We do not want to consider the bidomain source model from 2.2.1, because it needs data on the volume of the cardiac muscle. Moreover, as proved in [17, 18], throughout the myocardial volume an infinite dimensional distribution of transmembrane voltages can produce the same potential at the surface of the torso. Instead, we consider models that need only the heart and torso surfaces: the new EPI model from section 2.2.3, and the reference STD model from section 2.2.2. Hence, given a complete body surface data $z_T \in L^2(\Gamma_T)$, we consider three inverse problems as detailed in the next three sections.

The two first optimization problems being quadratic under a linear constraint, we use a Lagrangian formulation which allows to solve a linear system corresponding to the Lagrangian optimality conditions. We choose to state the optimization problems in the continuous setting and then discretize them.
### 2.3.1 STD source model with Tikhonov regularization (STD-Tk)

It consists of searching for $u_H \in H^1(\Gamma_H)$ that minimizes

$$J_{\text{std}}(u_H) = \frac{1}{2} \int_{\Gamma_T} |u(u_H) - z_T|^2 + \frac{\varepsilon}{2} \|\nabla u_H\|_{L^2(\Gamma_H)}^2, \quad (13)$$

where $u(u_H) \in H^1(\Gamma_H)/\mathbb{R}$ is the unique solution to the direct problem (8) with $\sigma_T \nabla u \cdot n_T = 0$ on $\Gamma_T$, and $u = u_H$ on the epicardium $\Gamma_H$. Let us consider a lifting operator $\mathcal{R}u_H \mapsto U_H$: such that $U_h \in H^1(\Omega_T)$, $U_H = u_H$ on $\Gamma_H$ and $\sigma_T \partial_n U_H = 0$ on $\Gamma_T$. We define $w = u - U_H$. With the definition of $U_H$ and the boundary condition for $u$ on $\Gamma_H$ we know that $w = 0$ on $\Gamma_H$, therefore $w$ belongs to $H_{0,\Gamma_H}^1 = \{v \in H^1(\Omega_T) \text{ s.t. } v|_{\Gamma_H} = 0\}$. The variational formulation of the problem (8) writes then:

$$\int_{\Omega_T} \sigma_T \nabla w \cdot \nabla \phi_0 = - \int_{\Omega_T} \sigma_T \nabla U_h \nabla \phi_0, \quad \forall \phi_0 \in H_{0,\Gamma_H}^1. \quad (14)$$

We thus define the Lagrangian function

$$\mathcal{L}(w, u_H, \lambda) = \frac{1}{2} \int_{\Omega_T} |w + U_H - z_T|^2 + \frac{\varepsilon}{2} \int_{\Gamma_H} \nabla u_H \cdot \nabla u_H + \int_{\Omega_T} \sigma_T \nabla w \cdot \nabla \lambda + \int_{\Omega_T} \sigma_T \nabla U_H \cdot \nabla \lambda, \quad (15)$$

with $\lambda \in H_{0,\Gamma_H}^1$. The optimality conditions can be written

$$0 = \langle \frac{\partial \mathcal{L}}{\partial w}, \phi_0 \rangle = \int_{\Gamma_T} (u - z_T) \phi_0 + \int_{\Omega_T} \sigma_T \nabla \lambda \cdot \nabla \phi_0, \quad \forall \phi_0 \in H_{0,\Gamma_H}^1, \quad (16a)$$

$$0 = \langle \frac{\partial \mathcal{L}}{\partial U_H}, \phi \rangle = \int_{\Gamma_T} (u - z_T) \phi + \varepsilon \int_{\Gamma_H} \nabla u_H \cdot \nabla \phi + \int_{\Omega_T} \sigma_T \nabla \lambda \cdot \nabla \phi, \quad \forall \phi \in H^1(\Omega_T) \quad (16b)$$

$$0 = \langle \frac{\partial \mathcal{L}}{\partial \lambda}, \phi_0 \rangle = \int_{\Omega_T} \sigma_T \nabla u \cdot \nabla \phi_0, \quad \forall \phi_0 \in H_{0,\Gamma_H}^1. \quad (16c)$$

These equations are then discretized using the same matrices and unknowns representation as described in section (2.2.2). Remark that some informations contained in equation (16a) are redundant with (16b). Indeed, the strong equations corresponding to (16a) and (16b) are respectively

$$\begin{align}
\sigma_T \Delta \lambda &= 0 && \text{in } \Omega_T, \nonumber \\
\sigma_T \partial_n \lambda + u &= z_T && \text{on } \Gamma_T, \nonumber \\
\lambda &= 0 && \text{on } \Gamma_H, \tag{17}
\end{align}$$

$$\begin{align}
\sigma_T \Delta \lambda &= 0 && \text{in } \Omega_T, \nonumber \\
\sigma_T \partial_n \lambda + u &= z_T && \text{on } \Gamma_T, \nonumber \\
\sigma_T \partial_n \lambda + \varepsilon \Delta u_H &= 0 && \text{on } \Gamma_H, \tag{18}
\end{align}$$

The only information from (17) not contained (18) is $\lambda = 0$ on $\Gamma_H$. Thus the whole equation (17) can be replaced by the condition $\lambda = 0$ on $\Gamma_H$ in the linear system to be solved.

Note that this source model is standard, as it is the most used for ECGi, but its finite-element discretization is more rarely employed than boundary elements or the method of fundamental solutions [2], as well as the first-order regularization. Moreover, the Lagrangian formulation of the optimization problem is also non-standard, as it amounts to optimize before discretize, while the most common approaches in ECGi rather use a discretize-then-optimize approach. However, these differences in the discretization process only produce small numerical variations with respect of the characteristic size of the meshes used for the inverse problem.

### 2.3.2 EPI source model with Tikhonov regularization (EPI-Tk)

It consists of searching for $\bar{v} \in H^1(\Gamma_H)$ that minimizes

$$J_{\text{epi-tk}}(\bar{v}) = \frac{1}{2} \int_{\Gamma_T} |u(\bar{v}) - z_T|^2 + \frac{\varepsilon}{2} \|\nabla \bar{v}\|_{L^2(\Gamma_H)}^2 + \frac{\varepsilon_{\text{inv}}}{2} \|\bar{v}\|_{L^2(\Gamma_H)}^2, \quad (19)$$

where $u(\bar{v})$ is obtained as the second component of the unique solution $(\bar{u}, u) \in V$ of (10) with $\sigma_T \nabla u \cdot n_T = 0$ on $\Gamma_T$. The term $\frac{\varepsilon_{\text{inv}}}{2} \|\bar{v}\|_{L^2(\Gamma_H)}^2$, with $\varepsilon_{\text{inv}} \ll \varepsilon$ is a penalty term that sets the
additive constant of $\bar{v}$ for practical computations. The actual regularization is provided by the term $\frac{\epsilon}{2}\|\nabla\bar{v}\|_{L^2(\Gamma_H)}^2$.

Note that, as compared to problem STD-Tk in section 2.3.1, the Tikhonov regularization applies to the epicardial transmembrane voltage $\bar{v}$ (instead of the electrical potential $u_H$) which has clearer electro-physiological meaning. Note in addition that the epicardial electrical potential is also reconstructed here as the trace on $\Gamma_H$ of $u$.

We thus define the Lagrangian function for $(\lambda_H, \lambda) \in (H^1(\Gamma_H) \times H^1(\Omega_T))$,

$$ \mathcal{L}(\bar{u}, u, v, \lambda_H, \lambda) = \frac{1}{2} \int_{\Omega_T} |u - z_T|^2 + \frac{\epsilon}{2} \int_{\Gamma_H} \nabla u \cdot \nabla v + \frac{\epsilon_{inv}}{2} \int_{\Gamma_H} |v|^2 + \int_{\Gamma_H} h(\bar{\sigma}_i + \bar{\sigma}_e) \nabla_s \bar{u} \nabla_s \lambda_H \\ + \int_{\Omega_T} \sigma_T \nabla u \cdot \nabla \lambda + \int_{\Gamma_H} \frac{\sigma^n}{\alpha h} (\bar{u} - u)(\lambda_H - \lambda) - \int_{\Gamma_H} h \bar{\sigma}_i \nabla_s v \nabla_s \lambda_H. \qquad (20) $$

The optimality conditions can be written

$$ 0 = \langle \frac{\partial \mathcal{L}}{\partial \bar{u}}, \phi_H \rangle = \int_{\Gamma_H} h(\bar{\sigma}_i + \bar{\sigma}_e) \nabla_s \phi_H \cdot \nabla_s \lambda_H + \int_{\Gamma_H} \frac{\sigma^n}{\alpha h} (\lambda_H - \lambda) \phi_H \quad \forall \phi_H \in H_{\Gamma_H}^1, \qquad (21a) $$

$$ 0 = \langle \frac{\partial \mathcal{L}}{\partial u}, \phi \rangle = \int_{\Omega_T} (u - z_T) \phi + \int_{\Omega_T} \sigma_T \nabla \phi \cdot \nabla \lambda - \int_{\Gamma_H} \frac{\sigma^n}{\alpha h} (\lambda_H - \lambda) \phi \quad \forall \phi \in H^1(\Omega_T), \qquad (21b) $$

$$ 0 = \langle \frac{\partial \mathcal{L}}{\partial v}, \phi_H \rangle = \varepsilon \int_{\Gamma_H} \nabla_s v \cdot \nabla_s \phi_H + \varepsilon_{inv} \int_{\Gamma_H} v \phi_H + \int_{\Gamma_H} h \bar{\sigma}_i \nabla_s \phi_H \cdot \nabla_s \lambda_H \quad \forall \phi_H \in H_{\Gamma_H}^1, \qquad (21c) $$

$$ 0 = \langle \frac{\partial \mathcal{L}}{\partial \lambda_H}, \phi_H \rangle = \int_{\Gamma_H} h(\bar{\sigma}_i + \bar{\sigma}_e) \nabla_s \bar{u} \cdot \nabla_s \phi_H + \int_{\Gamma_H} \frac{\sigma^n}{\alpha h} (\bar{u} - u) \phi_H \\ + \int_{\Gamma_H} h \bar{\sigma}_i \nabla_s \phi_H \cdot \nabla_s v \quad \forall \phi_H \in H_{\Gamma_H}^1, \qquad (21d) $$

$$ 0 = \langle \frac{\partial \mathcal{L}}{\partial \lambda}, \phi \rangle = - \int_{\Gamma_H} \frac{\sigma^n}{\alpha h} (\bar{u} - u) \phi + \int_{\Omega_T} \sigma_T \nabla u \cdot \nabla \phi \quad \forall \phi \in H^1(\Omega_T). \qquad (21e) $$

These equations are then discretized using the same matrices and unknowns representation as described in section 2.2.3.

In practice, the LU factorization of this system, once discretized, can be computed once and stored, enabling efficient resolution of the inverse problem over multiple time iterations. For an entire activation sequence, the reconstruction is fast (about 10 minutes for 140 time steps, including approximately 100 seconds for the LU factorization on a laptop with 32 GB RAM, using Python).

### 2.3.3 EPI source model with total variation regularization (EPI-TV)

It consists of searching for $\bar{v} \in H^1(\Gamma_H)$ that minimizes

$$ J_{\text{epi-tv}}(\bar{v}) = \int_{\Gamma_T} |u - z_T| + \varepsilon \|\nabla \bar{v}\|_{L^1(\Gamma_H)} + \frac{\varepsilon_{inv}}{2} \|\bar{v}\|_{L^2(\Gamma_H)}^2, \qquad (22) $$

where $u(\bar{v})$ is obtained as the second component of the unique solution $(\bar{u}, u) \in V$ of (10) with (10c) and $\sigma_T \nabla u \cdot n_T = 0$ on $\Gamma_T$. As compared to the EPI-Tk problem, the regularization term has changed, from $L^2$ on $\nabla \bar{v}$ to the total variation (TV) of $\bar{v}$. Indeed, it is well known that first-order Tikhonov regularization generate spurious oscillations, and excessively smooth solutions, while the transmembrane voltage is expected, for physiological reasons, to exhibit sharp wave fronts separating large regions with low variations. The STD source model does not provide the voltage $\bar{v}$, and in comparison, the new EPI source model allows TV regularization on the transmembrane voltage. Thus, it is expected to improve the quality of the reconstruction. In addition, the data fitting term has also been changed from an $L^2$ to an $L^1$ norm on the difference $u|_{\Gamma_T} - z_T$. This change is not motivated by the nature of the noise, taken Gaussian in our experiments, but rather by its numerical effect of automatically balancing the terms in Algorithm 1, leading to improved reconstruction quality.

In the (22) problem, the $J_{\text{epi-tv}}$ functional becomes non-linear and non-differentiable due to the presence of $L^1$ norms. The optimization problem is then solved using an iterative approach from [35], classically used in cardiology [37, 14] when regularizing on the $L^1$ norm of the u normal epicardial flux.
The "fixed point" algorithm [35] relies on the following approximation of the L¹ norm: for a small $\beta > 0$,

$$ \int_{\Gamma_H} |u| \approx \int_{\Gamma_H} \frac{u^2}{\sqrt{u^2 + \beta}} \approx \int_{\Gamma_H} \frac{u^2}{\sqrt{\tilde{u}^2 + \beta}}, $$

where $\tilde{u} \approx u$. The $\beta$ parameter is introduced to avoid the non-differentiability of the absolute value, while the replacement of $u$ in the denominator by a nearby function $\tilde{u}$ allows to linearize the approximate differential of the L¹ norm with respect to $u$. The iterative fixed point algorithm applied to the minimization problem (22) is described in Algorithm 1. At each iteration, we solve

**Algorithm 1 Fixed point method**

**Initialization:** $((\tilde{u}^0, u^0), v^0)$ solution of the Tikhonov epicardial model based ECGi problem EPI-Tk.

while not converged do

**Update:**

$$ ((\tilde{u}^{k+1}, u^{k+1}), v^{(k+1)}) := \underset{((\tilde{u}, u), v) \text{ satisfying (10)}}{\operatorname{arg min}} \frac{1}{2} \int_{\Gamma_T} \frac{(u - z_T)^2}{\sqrt{(u^k - z_T)^2 + \beta}} \\ + \frac{\varepsilon}{2} \int_{\Gamma_H} \frac{\nabla v \cdot \nabla v}{\sqrt{\|\nabla v^k\|_2^2 + \beta}} + \frac{\epsilon_{inv}}{2} \int_{\Gamma_H} v^2 $$

end while

**Output:** Solution $((\tilde{u}^{k+1}, u^{k+1}), v^{(k+1)})$

a quadratic subproblem, which entails the inversion of a Lagrangian optimality condition linear system.

# 3 Methodology for the evaluation of the forward EPI source model and inverse reconstructions

In this section we first present the general workflow for solving the direct and inverse problems and compute the activation maps. Then we describe the numerical methods and parameters used to produce both reference solutions and simulated datasets for the epicardial model and the ECGi inverse problem which results will be exposed in Sections 4 (forward validation of the EPI model) and 5 (results of the EPI-Tk, EPI-TV and STD ECGi problems).

## 3.1 General workflows

In the following, we assimilate the $\bar{v}$ source term of the EPI model to the trace of the transmembrane voltage on the epicardial surface $v|_{\Gamma_H}$.

### 3.1.1 Numerical validation of the EPI model

To numerically assess the EPI validity as a novel source model we compare, for several test case scenarios, solutions produced by the EPI model with solutions produced by the complete bidomain model, denoted as a reference $\cdot_{\text{ref}}$. For each test case, described in following subsections, the bidomain model produces both a transmembrane voltage $v_{\text{ref}}$ and an extracellular / extracardiac potential $u_{\text{ref}}$ temporal signals. The EPI model is a static application $\bar{v} \in H^1(\Gamma_H) \mapsto (\tilde{u}, u) \in H^1(\Gamma_H) \times H^1(\Omega_T)$, depending on a parameter set $(\alpha, p) \in ((0,1) \times (0,1))$. We recall that $p$ is such that $h(x) = p \times \text{wall thickness at } x$.

Denoting $\Pi_{\text{mesh}}$ a discrete interpolation operator between different meshes, we compare the trace on the epicardium $\Pi_{\Gamma_H} u_{\text{epi}}$ of the extracardiac potential $u_{\text{epi}}$ given by the EPI model $\Pi_{\text{mesh}}\Pi_{\Gamma_H} v_{\text{ref}} \mapsto (\tilde{u}, u_{\text{epi}})$ with the epicardial potential of reference $\Pi_{\text{mesh}}\Pi_{\Gamma_H} u_{\text{ref}}$, and so, for different parameter sets $(\alpha, p)$. In practice, projections (also called interpolations) $\Pi_{\text{mesh}}$ between meshes are computed using the RBFINTERPOLATOR operator from the SCIPY Python library. The validation workflow is summarized in Figure 3.
<figure><img src="image_5.png" /><figcaption>Figure 3: Numerical validation process for the forward EPI model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Workflow for solving the inverse ECGi problem. The reconstructed signals at each time step are then concatenated and post-processed to compute activation maps.</figcaption></figure>

### 3.1.2 Solving the static ECGi problems

We consider several simulated propagation scenarios, which potentials and voltages $u_{\text{ref}}$, $v_{\text{ref}}$ are computed by solving the anisotropic bidomain model on a torso-myocardium volume geometry. A second geometry, containing only the torso volume - and the epicardial surface considered as part of its boundary - is also considered. Simulated data for the ECGi problem are computed as noisy projections of the bidomain torso surface trace of the extracardiac potential on the torso-only mesh: $z_T = \Pi_{\text{mesh}}\Pi_{\Gamma_T}u_{\text{ref}} + \text{Gaussian noise}$. The Gaussian noise variance is chosen as 10% of the average amplitude in time of the original torso data. Then the three inverse problems EPI-Tk, STD-Tk and EPI-TV as defined in subsections 2.3.1, 2.3.2 and 2.3.3 are solved on the torso-only geometry for each time step of each propagation scenario data. Note that employing different meshes and geometries for data simulation and inverse resolution allow to avoid the so called "inverse crime". In addition, for EPI-based ECGi, the models employed for forward and inverse simulation are different and further avoid committing inverse crime. The methodology is summarized in Figure 4.

**Choice of regularization parameters for EPI-Tk, STD-Tk and EPI-TV.** For all three minimization problems, the regularization parameter $\varepsilon$ was selected using the L-curve method. For the EPI-Tk and STD-Tk formulations, the parameter obtained from the L-curve was increased by one order of magnitude in order to improve robustness to noise and reduce Tikhonov-induced spurious oscillations. For the EPI-TV problem and for the initialization phase of Algorithm 1, $\varepsilon$ was
<figure><img src="image_7.png" /><figcaption>Figure 5: L-curve associated to the inverse problem EPI-Tk (19). The optimal regularization parameter $\epsilon_{opt}$ gives the maximal curvature point of the presented curve. In our case, we have a range of optimal parameters such that $10^{-6} < \epsilon_{opt} < 5.10^{-5}$.</figcaption></figure>

chosen equal to the L-curve optimal parameter. This setup allows the algorithm to localize more quickly the activation front. An example of the L-curve obtained with the EPI-Tk problem is given Figure 5. The penalty factor $\epsilon_{inv}$ is set to $10^{-8}$ for both EPI-TV and EPI-Tk problems.

### 3.1.3 Activation maps computation

To reconstruct activation maps, inverse reconstructions obtained from data $z_T$ at successive time steps are concatenated. Then, post-processing algorithms are applied to the temporal reconstructions of $v$ and $u$. These post-processing steps are detailed below.

1. **(Pre-processing of the transmembrane potential)** At each time step, the potential $v$ can only be reconstructed up to an additive constant. This constant must therefore be adjusted to ensure a consistent temporal sequence. Moreover, due to the spatial regularization applied in the inverse problem, the amplitude of the reconstructions may vary over time, whereas true transmembrane voltages have fixed amplitude throughout the whole activation sequence. In particular, at the beginning (resp. end) of the activation sequence, the amplitude of reconstructed $v$ increases (resp. decreases) with data intensity, instead of reflecting an increasing (resp. decreasing) size of the activated (resp. resting) region of the heart.

To correct for this, two processing steps are applied to the temporal reconstructions of $v$. First, as in [29], the spatial median absolute deviation of $v$ at each time step is computed. The time of maximal median absolute deviation ($t_d$) indicates the midpoint of the depolarization sequence. Then, for $t < 2/3t_d$, the constant of $v$ is adjusted so that $\min v = 0$, while the amplitude of $v$ is left unchanged. This allows to preserve amplitude modulation induced by the regularization. For $t \ge 2/3t_d$, we apply a rescaling to enforce constant amplitude:

$$ \tilde{v}(t) = [\max_x(v(t_d)) - \min_x(v(t_d))] \frac{v(t) - \min_x(v(t))}{\max_x(v(t)) - \min_x(v(t))} $$

2. **(Temporal smoothing)** The temporal signals $u$, $v$, $\nabla u$, and $\nabla v$ are smoothed with a Gaussian filter of order 0 and standard deviation $\sigma = 10$ ms. Further investigation may be needed to determine an optimal filtering strategy, but was not the focus of this study.

3. **(Spatio-temporal derivative)** As in [30], the spatio-temporal derivative is defined as the product of the temporal derivative and the norm of the surface gradient of a signal. To compute the surface gradients of $u$ and $v$, we exploit the properties of $P^1$ Lagrange finite elements, which are affine within each triangle of the mesh. Due to Tikhonov regularization, gradients are naturally smooth. Temporal derivatives are computed with a centered Euler scheme applied to the time-smoothed values of $u$ and $v$. Activation times are then defined as the instant of maximal positive spatio-temporal derivative for $v$, and as the instant of maximal negative spatio-temporal derivative for $u$.
## 3.2 Meshes

### 3.2.1 Forward bidomain problem

<figure><img src="image_8.png" /><figcaption>Figure 6: Mesh used for bidomain simulations, with the myocardium mesh on the left and torso surface mesh on the right.</figcaption></figure>

All meshes were computed from an initial basis coarse mesh of a patient from Bordeaux University Hospital, describing the myocardium, blood cavities and torso geometry. This coarse mesh was refined and smoothed to produce the so called "bidomain mesh". The refinement was carried out using the Mmg software (see Figure 6). The bidomain mesh contains 17,523 nodes in the myocardium and 2,882 nodes on the torso surface.

### 3.2.2 EPI and STD problems

A torso-only volumic mesh derived from the same segmentation geometry was used for solving EPI and STD models. The boundaries of the torso mesh correspond to the epicardium on one side, artificially closed by valves, and to the external torso surface on the other (see Figure 7). This mesh contains 17,218 nodes on the epicardium ($\Delta x \approx 2$ mm), 16,220 nodes in the torso volume, and 1,475 nodes on the torso surface ($\Delta x \approx 30$ mm). Note that the valves used to artificially close

<figure><img src="image_9.png" /><figcaption>Figure 7: Mesh used to solve the equations of the epicardial (EPI) and standard (STD) models. The torso surface is shown on the right, and the epicardium (red) artificially closed by valves (orange) on the left.</figcaption></figure>

the epicardial mesh are treated numerically as part of the torso domain. Homogeneous Neumann conditions are therefore imposed on fluxes across the valves boundaries in the FEM variational formulation of the EPI model.

## 3.3 Bidomain parameters

### 3.3.1 Ionic model and electrical stimulation

The parameters of the Mitchell-Shaeffer model are set as in [24]. The stimulus current (pacing) $s(x, t)$ is used to locally stimulate the heart and initiate the propagation of an action potential. We model $s(x, t)$ as the product of a spatial smoothed rectangular pulse - the 0-1 transition is smoothed by a cubic function - and of a temporal rectangular pulse.
<figure><img src="image_10.png" /><figcaption>Figure 8: Fiber directions included in the anisotropic bidomain simulations. The function $\phi_{\text{epi}}$ (color scale) indicates the heart depth: blue arrows represent endocardial fibers while red arrows represent epicardial fibers. The fiber directions were computed following the rule from [1].</figcaption></figure>

### 3.3.2 Cardiac fibers and conductivities

For anisotropic propagation scenarios, fiber directions were generated using the rule-based method of Bayer et al. [1]. An illustration of computed fiber directions is shown Figure 8. In the case of anisotropic conductivity tensors, the eigenvalues were taken from the work of Potse [25]. Denoting $R$ the transition matrix to the fiber direction basis, the heart conductivity tensors write

$$ \sigma_i = R \begin{pmatrix} 3 & 0 & 0 \\ 0 & 0.3 & 0 \\ 0 & 0 & 0.3 \end{pmatrix} R^T, \quad \sigma_e = R \begin{pmatrix} 3 & 0 & 0 \\ 0 & 1.2 & 0 \\ 0 & 0 & 1.2 \end{pmatrix} R^T. $$

For isotropic cases, the conductivity tensors in the heart were chosen as $\sigma_i = I_3$ and $\sigma_e = 3I_3$. In both isotropic and anisotropic cases, the torso conductivity was set to $\sigma_T = 2$.

## 3.4 ECGi source models parameters

### 3.4.1 Conductivities

The conductivity tensor $\sigma_T$ appearing in the STD model is chosen identical to the bidomain torso conductivity tensor.

For isotropic propagation scenarios, the conductivities of the EPI model are chosen identical to the conductivities of the bidomain model. For anisotropic scenarios, specifically for ECGi reconstructions, we take $\bar{\sigma}_i = 2I_2$, $\bar{\sigma}_e = 3I_2$, $\sigma_T = 2I_3$, and $\sigma^n = 1.5$. These conductivity values are chosen arbitrarily, while remaining within ranges consistent with the eigenvalue spectrum of the bidomain conductivity tensors.

### 3.4.2 Computation of the EPI thickness parameter

The thickness parameter $h$ used in the EPI novel source model 2.2.3 is defined as a fixed fraction $p$ of the ventricular wall thickness $h_{\text{full}}$: $h = p h_{\text{full}}$. Since the ventricular wall thickness varies spatially (the right ventricle is approximately half as thick as the left ventricle), $h_{\text{full}}$ must be defined as a function of $x$, also discretized on the $P^1$ Lagrange basis. The full ventricular wall thickness $h_{\text{full}}$ is computed as a distance from the epicardium to the endocardium on a fine myocardial mesh, using the *Heat method* [8]. The complete algorithm is detailed in [20].

On the ventricles, $h_{\text{full}}$ closely matches the true ventricular wall thickness. However, above the septum, the wall thickness is not well-defined, and $h_{\text{full}}$, constructed as a distance to the endocardium, behaves like a local average thickness over a small surface.

### 3.4.3 EPI parameters ($p, \alpha$)

In the EPI model validation workflow, a wide range of parameter sets $(p, \alpha) \in ((0, 1) \times (0, 1))$ is taken to both evaluate performances of the EPI model and calibrate $p$ and $\alpha$. For EPI-based ECGi
<figure><img src="image_11.png" /><figcaption>Figure 9: Pacing sites, from left to right: left-ventricle (LV), right-ventricle (RV), Apex (APEX), Right-ventricular-Outflow-Tract (RVOT)</figcaption></figure>

reconstructions, an optimal parameter set $(p, \alpha)$ is chosen accordingly to the calibration study, see section 4.

## 3.5 Table of data scenarios

The simulated data scenarios are summarized in the Table 1. A scenario is defined by the choice of one or two stimulation sites, see Figure 9, with location indicated in the table, and the use of isotropic or anisotropic cardiac conductivity tensors as defined in previous subsection 3.3. Isotropic (resp. anisotropic) data scenarios will be denoted with their names and the indication "iso" (resp "aniso"). For scenarios with two pacing sites (RV+APEX and RVOT+APEX), the delay between the two pacing times is specified in Table 1.

<table><caption>Table 1: Names and pacing characteristics of the data scenarios. All protocols are simulated with two propagation types: isotropic and anisotropic.</caption><thead><tr><th>Name</th><th>Pacing site 1</th><th>Pacing site 2</th><th>Delay (t2 - t1)</th></tr></thead><tbody><tr><td>LV</td><td>LV</td><td>-</td><td>-</td></tr><tr><td>RV</td><td>RV</td><td>-</td><td>-</td></tr><tr><td>APEX</td><td>APEX</td><td>-</td><td>-</td></tr><tr><td>RVOT</td><td>RVOT</td><td>-</td><td>-</td></tr><tr><td>RV+APEX</td><td>APEX</td><td>RV</td><td>20 ms</td></tr><tr><td>RVOT+APEX</td><td>APEX</td><td>RVOT</td><td>20 ms</td></tr></tbody></table>

For the validation of the EPI forward model, data generated from both isotropic ("iso") and anisotropic ("ani") scenarios are considered. For convenience, we will hereafter refer to data generated under isotropic or anisotropic conditions simply as isotropic or anisotropic data, respectively. For the ECGi reconstructions, only data generated from anisotropic scenarios are used to evaluate the performance of the EPI-Tk, EPI-TV, and STD-Tk optimization problems.

## 3.6 Evaluation metrics

For simplicity of notation, we refer to the interpolated reference signals as $u_{\text{ref}}$, $v_{\text{ref}}$ instead of $\Pi_{\text{mesh}}u_{\text{ref}}$ and $\Pi_{\text{mesh}}v_{\text{ref}}$ in the following.

The metrics used in this paper are derived from two classical metrics: the correlation coefficient and the $L^2$ relative error. For a size $N_{\text{nodes}}$ vector $\mathbf{X}$ defined on a mesh, we denote $\mathbf{X}_i, i \in [1, N_{\text{nodes}}]$ the value of $\mathbf{X}$ at the $i^{\text{th}}$ node of the mesh, and $\bar{\mathbf{X}}$ its average value over all nodes. The correlation coefficient between $\mathbf{X}$ and a reference $\mathbf{X}_{\text{ref}}$ writes

$$CC = \frac{\sum_{i=1}^{N_{\text{nodes}}} (\mathbf{X}_i - \bar{\mathbf{X}})(\mathbf{X}_i - \bar{\mathbf{X}}_{\text{ref}})}{\sqrt{\sum_{i=1}^{N_{\text{nodes}}} (\mathbf{X}_i - \bar{\mathbf{X}})^2} \sqrt{\sum_{i=1}^{N_{\text{nodes}}} (\mathbf{X}_{\text{ref},i} - \bar{\mathbf{X}}_{\text{ref}})^2}} \quad (23)$$

The $L^2$ norm relative error of $\mathbf{X}$ with respect to a reference $\mathbf{X}_{\text{ref}}$ writes

$$\frac{\|\mathbf{X} - \mathbf{X}_{\text{ref}}\|_{L^2}}{\|\mathbf{X}_{\text{ref}}\|_{L^2}} = \frac{\sqrt{(\mathbf{X} - \mathbf{X}_{\text{ref}})^T M (\mathbf{X} - \mathbf{X}_{\text{ref}})}}{\sqrt{\mathbf{X}_{\text{ref}}^T M \mathbf{X}_{\text{ref}}}} \quad (24)$$
where M is the P¹ Lagrange finite element mass matrix on the considered mesh.

### 3.6.1 EPI model validation

In the validation workflow described in Figure 3, three types of measurements are employed to compare the EPI model potentials with their reference counterparts:

* **Temporal cross-correlation $\overline{CC}$**: For each time step, we compute the correlation coefficient $CC(t = t_k)$ between the epicardial model potential $u_{\text{epi}}$ on the epicardium and the reference potential $u_{\text{ref}}$ interpolated on the epicardial mesh. Then these correlation coefficients are averaged over all considered time steps $t_k$, $k \in [1, N_{\text{it}}]$:

$$ \overline{CC} := \frac{1}{N_{\text{it}}} \sum_{k=1}^{N_{\text{it}}} CC(t = t_k). \qquad (25) $$

* **Normalized $L^2$ relative error $||u_{\text{epi}}^{\text{scaled}} - u_{\text{ref}}^{\text{scaled}}||/||u_{\text{ref}}^{\text{scaled}}||$**: Across all time steps, the mean values of $u_{\text{epi}}$ and $u_{\text{ref}}$ are shifted to zero, and both signals are rescaled to lie within [0, 1] over the full time interval. Explicitly, we apply

$$ u^{\text{scaled}}(t, x) = \frac{u(t, x) - \min_{t,x} u(t, x)}{\max_{t,x} u(t, x) - \min_{t,x} u(t, x)}. \qquad (26) $$

At each time step, we compute the relative $L^2$ error between $u_{\text{epi}}^{\text{scaled}}$ and $u_{\text{ref}}^{\text{scaled}}$, and finally average these errors over time.

* **Averaged $L^2$ relative error $||u_{\text{epi}} - u_{\text{ref}}||/||u_{\text{ref}}||$**: Across all time steps, the mean values of $u_{\text{epi}}$ and $u_{\text{ref}}$ are again shifted to zero, then the relative $L^2$ error between $u_{\text{epi}}$ and $u_{\text{ref}}$ are computed for each time step and averaged over time.

### 3.6.2 ECGi reconstructions

**Quantitative assessment** The temporal cross correlation $\overline{CC}$ metric is also used for the assessment of reconstructed transmembrane voltage $v(x, t)$ and epicardial extracellular potential $u(x, t)$.

For activation maps, three metrics are considered: the relative $L^2$ error, the correlation coefficient $CC$, and the Slowness Coefficient $SC$. The Slowness Coefficient introduced in [30] corresponds to the correlation coefficient of the surface gradient of the activation map. This metric is particularly useful as it is sensitive to block lines, unlike the standard correlation coefficient.

**Qualitative assessment** We also use qualitative visual assessment for both activation maps and reconstructed $v$ and $u$. Specifically we focus on

* the shape of the activation front at some fixed time steps;

* local electrograms and their transmembrane voltage equivalent;

* activation maps features such as artificial lines of block or earliest activation sites.

## 4 Forward numerical assessment of the EPI model

In this section, we present the results of the validation workflow detailed in Section 3 and summarized in Figure 3.

### 4.1 Calibration

In this subsection, only the isotropic scenarios from Table 1 are studied, with the purpose of assessing the sensitivity of the epicardial model to variations in $p$ and $\alpha$ parameters, and find values that minimize the discrepancy with bidomain model solutions. In such isotropic setting, the effective intracellular, effective extracellular, and torso conductivities from the epicardial model can be chosen identical to those of the reference bidomain simulation (see Section 3). The sensitivity of the model to the value of $(\alpha, p)$ is illustrated in Figure 10, showing the EPI model solution $u_{\text{epi}}$ associated to various sets $(\alpha, p)$ for a given time $t$. In this example, the overall spatial pattern of $u_{\text{epi}}$ remains close to that of the reference solution $u_{\text{ref}}$, but the amplitude of the potential varies significantly with
<figure><img src="image_12.png" /><figcaption>Figure 10: Examples of epicardial potentials $u_{\text{epi}}$ obtained with the epicardial model for different parameter sets $(\alpha, p)$, compared with the reference potential $u_{\text{ref}} = u_{\text{bid}}$ on the epicardium, for scenario RVOT+APEX iso.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Errors between the epicardial model potential $u_{\text{epi}}$ on the epicardium and the reference potential $u_{\text{ref}}$ in the $(p, \alpha)$ plane, for scenario LV iso.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 12: Errors between the epicardial (EPI) model potential $u_{\text{epi}}$ on the epicardium and the reference potential $u_{\text{ref}}$ in the $(p, \alpha)$ parameter plane, for scenario RVOT+APEX iso.</figcaption></figure>
the EPI model parameters. This observation is confirmed by the behavior of the three validation error norms introduced in Section 3.6 across the (p, α) parameter plane, for two scenarios: LV iso, and RVOT+APEX iso (Figures 11 and 12). For both scenarios, a consistent pattern emerges: the temporal correlation remains high (> 0.9) regardless of the parameter pair, and the normalized relative L² error (||u<sub>epi</sub><sup>scaled</sup> - u<sub>ref</sub><sup>scaled</sup>||/||u<sub>ref</sub><sup>scaled</sup>||) stays of the order of 10% throughout the whole parameter space. Note that while this error may appear to be large, it is in fact an acceptable ratio in the context of ECGi inverse problem. Indeed, in a clinical setting, other sources of modeling error, such as neglecting the lungs or cardiac fibers, can introduce much larger discrepancies. Then, the coupling of a high correlation combined with a low normalized L² error indicates that the epicardial model captures the spatial distribution of the bidomain solutions. By contrast, the non-normalized L² errors are larger, reflecting mainly amplitude discrepancies. With the chosen scaling (26) on epicardial potentials, the low normalized L² relative errors also indicates that this amplitude error is predominantly constant over time.

For the three inverse problem formulations presented in this article, an inaccurate representation of extracardiac potential amplitude is not critical. Because the applied regularizations (Tk, TV) do not constrain the amplitude of the reconstructed transmembrane potential, the amplitude of v can naturally adjust to fit the data. Consequently, the absolute amplitude of the solutions is of limited importance, provided it is consistent over time, as shown here, and the spatial distribution and activation front are accurately captured.

Nevertheless, across all error norms, an optimal parameter set emerges around (p, α) = (0.7, 0.1).
This parameter choice will be adopted for solving the inverse problems constrained by the epicardial
(EPI) source model.

## 4.2 Validation for anisotropic propagation scenarios

In order to further assess the use of the EPI model for ECGi, we address anisotropic propagation
test cases. Since the EPI model assumes isotropic surface conductivities, it introduces an additional
modeling error. This setting thus aims at validating that the EPI model can still be employed under
more realistic conditions where electrical conductivities and fiber orientation are not precisely known.
We emphasize in particular on amplitude-independent errors (CC and normalized L² relative errors),
in order to assess the model’s ability to reproduce the correct shape of the extracardiac potential
despite modeling inaccuracies.

<figure><img src="image_15.png" /><figcaption>Figure 13: Errors between the epicardial model potential $u_{\text{epi}}$ on the epicardium and the reference potential $u_{\text{ref}}$ in the $(p, \alpha)$ plane, for scenario LV aniso.</figcaption></figure>

For the two cases shown in Figures 13 and 14, the EPI model remains able to accurately reproduce
the shape of the extracardiac potential on the epicardium. Indeed, temporal cross-correlation values
are high across the entire parameter plane, and the errors $||u_{\text{epi}}^{\text{scaled}} - u_{\text{ref}}^{\text{scaled}}||/||u_{\text{ref}}^{\text{scaled}}||$ remain of the
order of 10%, with minima still attained around $(p, \alpha) = (0.7, 0.1)$.

**Remark 1** The errors reported here, especially those depending on the amplitude of the solution, should be interpreted with caution, as they are sensitive to the surface conductivity values chosen in the epicardial model. More importantly, the results indicate that, irrespective of the amplitude and the chosen parameter values, the epicardial model is able to reproduce the extracardiac potential behavior observed with the bidomain model.
<figure><img src="image_16.png" /><figcaption>Figure 14: Errors between the epicardial model potential $u_{epi}$ on the epicardium and the reference potential $u_{ref}$ in the $(p, \alpha)$ plane, for scenario RVOT+APEX aniso.</figcaption></figure>

# 5 Results for the inverse ECGi problem

This section is devoted to the evaluation of reconstructions performed with the three inverse problems STD-Tk, EPI-Tk and EPI-TV, following the workflow detailed in section 3, and summarized in Figure 4. We recall that all data scenarios used in this section are anisotropic.

## 5.1 Activation maps reconstructions

### 5.1.1 Quantitative assessment

<table><caption>Table 2: Errors on the activation maps reconstructed with solutions of the EPI-Tk problem ($\partial_{st}v_{Tk}$ and $\partial_{st}u_{Tk}$), of the EPI-TV problem ($\partial_{st}v_{TV}$ and $\partial_{st}u_{TV}$) and of the STD-Tk problem ($\partial_{st}u_{STD}$). The best values in each column (=test case) are highlighted.</caption><thead><tr><th>Method</th><th>Error</th><th>LV</th><th>RV</th><th>APEX</th><th>RVOT</th><th>RV+APEX</th><th>RVOT+APEX</th></tr></thead><tbody><tr><td rowspan="3">∂<sub>st</sub>v<sub>Tk</sub></td><td>L<sup>2</sup></td><td>0.08</td><td>0.08</td><td>0.09</td><td>0.10</td><td>0.07</td><td>0.16</td></tr><tr><td>CC</td><td>0.98</td><td>0.98</td><td>0.98</td><td>0.98</td><td>0.99</td><td>0.94</td></tr><tr><td>SC</td><td>0.96</td><td>0.96</td><td>0.98</td><td>0.95</td><td>0.98</td><td>0.94</td></tr><tr><td rowspan="3">∂<sub>st</sub>u<sub>Tk</sub></td><td>L<sup>2</sup></td><td>0.08</td><td>0.10</td><td>0.06</td><td>0.09</td><td>0.10</td><td>0.12</td></tr><tr><td>CC</td><td>0.98</td><td>0.98</td><td>0.99</td><td>0.98</td><td>0.97</td><td>0.96</td></tr><tr><td>SC</td><td>0.93</td><td>0.96</td><td>0.99</td><td>0.89</td><td>0.96</td><td>0.93</td></tr><tr><td rowspan="3">∂<sub>st</sub>v<sub>TV</sub></td><td>L<sup>2</sup></td><td>0.07</td><td>0.09</td><td>0.08</td><td>0.12</td><td>0.10</td><td>0.15</td></tr><tr><td>CC</td><td>0.99</td><td>0.98</td><td>0.99</td><td>0.96</td><td>0.98</td><td>0.95</td></tr><tr><td>SC</td><td>0.98</td><td>0.95</td><td>0.99</td><td>0.93</td><td>0.98</td><td>0.93</td></tr><tr><td rowspan="3">∂<sub>st</sub>u<sub>TV</sub></td><td>L<sup>2</sup></td><td>0.07</td><td>0.07</td><td>0.05</td><td>0.08</td><td>0.08</td><td>0.09</td></tr><tr><td>CC</td><td>0.99</td><td>0.99</td><td>0.99</td><td>0.98</td><td>0.99</td><td>0.98</td></tr><tr><td>SC</td><td>0.97</td><td>0.98</td><td>0.99</td><td>0.95</td><td>0.98</td><td>0.97</td></tr><tr><td rowspan="3">∂<sub>st</sub>u<sub>std</sub></td><td>L<sup>2</sup></td><td>0.06</td><td>0.09</td><td>0.06</td><td>0.08</td><td>0.08</td><td>0.11</td></tr><tr><td>CC</td><td>0.99</td><td>0.98</td><td>0.99</td><td>0.99</td><td>0.99</td><td>0.97</td></tr><tr><td>SC</td><td>0.97</td><td>0.97</td><td>0.99</td><td>0.96</td><td>0.98</td><td>0.95</td></tr></tbody></table>

We first evaluate the reconstructed activation maps for the three problems EPI-Tk, EPI-TV and STD-Tk through the error metrics defined in section 3. The errors for all anisotropic test cases are presented in Table 2. Considering the $L^2$, CC and SC norms, all methods appear comparable in terms of accuracy. Most of the time, the errors for all methods remain within a 0.01 interval, which is not significant given the level of noise. Consequently, on a quantitative basis, no predominant methodology can be extracted from these experiments, while, because of the addition of a physiological regularization, it would have been expected for the EPI-TV method to improve the accuracy of reconstructed activation maps compared to both EPI-Tk and STD-Tk.

A more qualitative analysis is then required to compare more deeply the three methods.
5.1.2 Qualitative assessement

**A- Tikhonov-regularized solutions** The qualitative assessment of reconstructed activation maps is based on two features: the presence of artificial lines of block, a long standing issue in ECGi [30], and the visual identification of earliest activation sites. Figure 15 presents activation maps

<figure><img src="image_17.png" /><figcaption>Figure 15: Activation maps (ms) reconstructed using the EPI-Tk problem (with *u* and *v*) and the STD-Tk problem. The stimulation sites for the four examples are shown with a red mark.</figcaption></figure>

reconstructed with the EPI-Tk and STD-Tk optimization problems, with a particular emphasis on stimulation sites (or earliest activation sites). It can be observed that both source models EPI and STD allow to identify visually the stimulation sites.

<figure><img src="image_18.png" /><figcaption>Figure 16: Activation maps (ms) reconstructed using the EPI-Tk problem (with *u* and *v*) and the STD-Tk problem. Some artificial lines of block (LOBs) are circled in red.</figcaption></figure>

Figure 16 displays the same activation maps, but with a focus on the presence of artificial lines of block (LOBs). Again, both STD-Tk and EPI-Tk optimization problems generate artificial LOBs, in similar or different heart locations.
Echoing the quantitative section, we observe that the sole error norms CC, L² and SC do not offer a general overview of reconstructed activation maps accuracy. Taking specific case RVOT+APEX, error norms suggest the superiority of STD-Tk reconstructions over EPI-Tk. However, qualitatively, both reconstructions are of equally poor quality, and present a strong, very visible artificial LOB. This suggests the need for future developments in quantification tools for activation maps assessment.

**Discussion: sole effect of using the EPI model** For both STD-Tk and EPI-Tk optimization problems, we observe a comparable accuracy of the reconstructed activation maps. On the one hand, this suggests that the sole use of the epicardial source model in the inverse problem does not yield a significant improvement in activation sequence recovery. On the other hand, these results imply that the introduction of modeling errors inherent to the epicardial model does not degrade the accuracy of the inverse reconstructions, especially in anisotropic scenarios. We emphasize that, except from the geometry disparities between forward and inverse simulations, the STD model does not include modeling error. Then, with the same amount of knowledge introduced in the models (through conductivities a priori and regularization), and despite modeling errors, the EPI model based reconstructions are of comparable accuracy to those of the Laplace equation based.

Additionally, the EPI model offers certain advantages over the classical STD one. It provides simultaneous access to both the transmembrane voltage and the extracellular potential on the epi-cardium, and allows the introduction of heterogeneity or anisotropy in the epicardial conductivity tensors, while being substitutable to the STD model in clinical devices. It opens new possibilities for introducing constraints in the ECGi problem, either through advanced regularization strategies or by incorporating prior knowledge on conductivity on the heart surface - incorporating scar knowledge for instance.

**B- Addition of a TV regularization** The EPI-TV optimization problem is precisely one example of such information incorporation, through the addition of a physiological total variation (TV) regularization on the transmembrane voltage. Figure 17 presents examples of reconstructed activation maps derived from solutions *u* and *v* of the EPI-TV problem. The red circles indicate the localization of former EPI-Tk artificial LOBs, that disappeared with the TV regularization.

<figure><img src="image_19.png" /><figcaption>Figure 17: Activation maps reconstructed from *v* and *u* solutions of the EPI-TV problem. The red circles indicate the localization of former EPI-Tk artificial LOBs, that disappeared with the TV regularization.</figcaption></figure>

Qualitatively, these are smoother than the EPI-Tk and STD-Tk recovered activation maps (AT maps), shown in Figure 16, as some artificial LOBs have been eliminated by the TV regularization. Examples of such smoothing are circled in red in Figure 17.
## Conclusion on AT maps recovery

While quantitative metrics do not clearly identify one methodology as superior, qualitative assessment indicates that the EPI-TV method outperforms the EPI-Tk and STD-Tk approaches, both of comparable accuracy. Specifically, the EPI-TV method eliminates some artificial lines of block (LOBs) that might otherwise falsely suggest the presence of pathological conduction abnormalities.

## 5.2 Signals reconstruction

This section is devoted on evaluating reconstructions of transmembrane voltage *v* and extracellular potential *u* signals with the EPI-Tk, STD-Tk and EPI-TV methods. Table 3 reports average correla-

<table><caption>Table 3: Average correlation coefficients over time of $v_{TV}$, $u_{TV}$, $v_{Tk}$, $u_{Tk}$, and $u_{STD}$, obtained by solving the EPI-TV, EPI-Tk and STD-Tk problems respectively, with respect to the references $v_{ref}$ and $u_{ref}$. The highest values for *v* and *u* in each row are highlighted.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Average temporal correlation coefficient CC</th></tr><tr><th>v<sub>TV</sub></th><th>v<sub>Tk</sub></th><th>v<sub>TV</sub></th><th>v<sub>Tk</sub></th><th>u<sub>STD</sub></th></tr></thead><tbody><tr><th>LV</th><td>0.77</td><td>0.74</td><td>0.62</td><td>0.51</td><td>0.52</td></tr><tr><th>RV</th><td>0.79</td><td>0.78</td><td>0.60</td><td>0.49</td><td>0.51</td></tr><tr><th>APEX</th><td>0.80</td><td>0.77</td><td>0.61</td><td>0.49</td><td>0.48</td></tr><tr><th>RVOT</th><td>0.66</td><td>0.67</td><td>0.54</td><td>0.43</td><td>0.43</td></tr><tr><th>RV+APEX</th><td>0.75</td><td>0.75</td><td>0.52</td><td>0.43</td><td>0.43</td></tr><tr><th>RVOT+APEX</th><td>0.71</td><td>0.70</td><td>0.49</td><td>0.40</td><td>0.40</td></tr></tbody></table>

tion coefficients over all time steps between reconstructions $v_{TV}$, $v_{Tk}$, $u_{TV}$, $u_{Tk}$ and $u_{STD}$ and their respective references, for the test activation sequences described in Section 3. Quantitatively, reconstructions $v_{TV}$ and $u_{TV}$ outperform $v_{Tk}$, $u_{Tk}$ and $u_{STD}$, as expected, with the most significative improvements on the extracellular potential *u*. Additionally, across all tested activation sequences, the reconstructions $u_{Tk}$ and $u_{STD}$ admit comparable correlation coefficients with the reference, in agreements with activation maps reconstruction results.

### 5.2.1 Signals space distribution

To visualize the TV improvement, Figure 18 illustrates examples of transmembrane voltages and extracellular potentials reconstructed at fixed time instants with problems EPI-TV ($v_{TV}$, $u_{TV}$) and EPI-Tk ($v_{Tk}$, $u_{Tk}$) initialization step of Algorithm 1. Qualitatively, away from the activation front, the reconstructions $u_{TV}$ and $v_{TV}$ exhibit reduced spurious oscillations compared to $u_{Tk}$ and $v_{Tk}$. In particular, the transmembrane voltage $v_{TV}$ displays more clearly delimited plateau regions, while the extracellular potential $u_{TV}$ shows physiologically smooth spatial variations within these regions. Simultaneously, the activation front itself is sharper.

### 5.2.2 Local electrograms (EGMs)

We plot local electrograms and preprocessed transmembrane voltage time course for three randomly chosen nodes on the epicardium, $n_1$, $n_2$, $n_3$, for scenario RV+APEX (Figure 19). Concerning *v* temporal courses, the EPI-TV produces meaningful action potentials, more consistent with the reference than EPI-Tk-based action potentials. This observation extends to local electrograms (EGMs). Specifically, in subfigure 19f, the $u_{EPI-TV}$ time courses follows almost perfectly the reference curve, while both EPI-Tk and STD-Tk electrograms do not reach the reference extremums. Additionally, in both subfigures 19d and 19f, the EPI-Tk EGMs are visually more accurate than the STD-Tk EGMs.

## Conclusion on signals recovery

In the case of signals recovery, both quantitative and qualitative assessment agree on the superiority of the EPI-TV method against both the EPI-Tk and STD methods. For local ECGs, it is also worth noting that the extracellular potential time courses obtained with EPI-Tk appear to follow the reference more closely than those obtained with the STD method.
<figure><img src="image_20.png" /><figcaption>Figure 18: Comparison at four distinct time instants ((a), (b), (c), (d)) of the reconstructed trans-membrane voltages *v* and extracellular potentials *u* obtained with optimization problems EPI-Tk and EPI-TV, against the reference solutions *v*<sub>ref</sub> and *u*<sub>ref</sub>.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(a) In time $v_{\text{ref}}$, $v_{\text{EPI-TV}}$ and $v_{\text{EPI-Tk}}$ at node $n_1$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(b) In time $v_{\text{ref}}$, $v_{\text{EPI-TV}}$ and $v_{\text{EPI-Tk}}$ at node $n_2$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c) In time $v_{\text{ref}}$, $v_{\text{EPI-TV}}$ and $v_{\text{EPI-Tk}}$ at node $n_3$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(d) Electrocardiograms for $u_{\text{ref}}$, $u_{\text{EPI-TV}}$, $u_{\text{EPI-Tk}}$ and $u_{\text{STD-Tk}}$ at node $n_1$.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(e) Electrocardiograms for $u_{\text{ref}}$, $u_{\text{EPI-TV}}$, $u_{\text{EPI-Tk}}$ and $u_{\text{STD-Tk}}$ at node $n_2$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(f) Electrocardiograms for $u_{\text{ref}}$, $u_{\text{EPI-TV}}$, $u_{\text{EPI-Tk}}$ and $u_{\text{STD-Tk}}$ at node $n_3$.</figcaption></figure>

<caption>Figure 19: In time reconstructions of *v* (preprocessed as detailed in section 3.1.3) and *u* on three different fixed points of the heart (from left to right).</caption>
**Remark 2 (Oscillations)** As a general comment, artifactual oscillations are a consequence of both Gaussian noise on the data, and Tikhonov regularization for Tk and STD signals. Indeed, Tikhonov regularization is well known to generate spurious unphysiological oscillations. When dealing with inverse problems, the regularization is always set as a compromise between obtaining good mathematical properties for the optimization problem, and not affecting too much the solution of this inverse problem. With unphysiological Tikhonov regularization, this compromise further includes dealing between seeing the very steep activation front in u and v and avoiding spurious oscillations. In the case of TV regularization on the transmembrane voltage v, the latter does not exist anymore as the regularization is really more physiological.

**Remark 3 (Practical usage of EPI-TV)** Our results indicate that solving the inverse problem with total variation regularization on the transmembrane voltage v enhances both the reconstruction of electrical potentials over the activation sequence and the quality of activation maps, which exhibit fewer artificial discontinuities. Nevertheless, solving the inverse problem with $L^1$ type norms is computationally more intensive than with Hilbert norms. The choice to adopt this approach should therefore depend on the intended objective. For example, activation maps reconstructed from the EPI-TV problem are smoother, but remain of the same order of accuracy as those from Tikhonov-based formulations. By contrast, when the objective is to directly reconstruct cardiac electrical signals, the total variation regularization on v, allowed by the use of the epicardial (EPI) model, provides a substantial and meaningful advantage.

# 6 Discussion

The epicardial (EPI) model constitutes an alternative source model for the ECGi inverse problem, which can be used in place of existing approaches. Its main advantage lies in coupling the two variables, *v* and *u*, on the epicardial surface with the cardiac volume, whereas classical models provide access to only one of these electrical signals. With a single source model, it is therefore possible to combine a wider range of regularization strategies or activation map computation methods, thereby extracting more information from the same data. Moreover, the epicardial model admits surface conductivities on the epicardium, which could be exploited to represent fibrotic or scarred regions. In that sense, the epicardial model provides a framework analogous to the bidomain electrostatic source model, while requiring only the torso geometry. It thus retains the simplicity necessary for practical use, while offering more modeling flexibility than, for instance, the Laplace-based formulation (STD model).

In this study based on simulated data, we have demonstrated that, under standard regularization schemes for the inverse problem, and no prior specific information, the use of the epicardial model within the ECGi framework produces activation map reconstructions of comparable accuracy to those obtained from solutions of the Laplace-constrained torso problem. This finding indicates that the modeling errors inherent to the epicardial formulation do not compromise the quality of the reconstructions.

As an example of application of the EPI model, we studied the effects of coupling the EPI model with a total variation regularization on the transmembrane voltage in ECGi reconstructions. Due to the very specific spatial shape of the transmembrane voltage, the total variation regularization on the transmembrane voltage is more physiological than applied on the extracellular potential. The results, based on simulated data, demonstrated a superiority of the EPI-TV method on signal reconstructions and electrograms of both *u* and *v*, and a reduction in artificial lines of block in the reconstructed activation maps.

Further investigations are nonetheless required to identify specific scenarios in which the epicardial model may provide more significant advantages over classical inverse problem formulations. Several research directions appear promising in this regard: for instance, using the epicardial model to validate or refute the presence of lines of block observed in activation maps, integrating prior knowledge about the presence of scars within the myocardium, or enforcing pathological patterns for the transmembrane voltage *v*. In addition, future work will focus on validating the use of the EPI model based ECGi in a more realistic clinical or experimental setting.

# Acknowledgments

Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux
INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux and CNRS (and ANR in accordance to the programme d'investissements d'Avenir (see http://www.plafrim.fr/)).

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] Jason D Bayer et al. "A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models". In: *Annals of biomedical engineering* 40 (2012), pp. 2243-2254.

[2] Laura Bear et al. "Introduction to noninvasive cardiac mapping". In: *Cardiac electrophysiology clinics* 7.1 (2015), pp. 1-16.

[3] Laura R Bear et al. "Advantages and pitfalls of noninvasive electrocardiographic imaging". In: *Journal of electrocardiology* 57 (2019), S15-s20.

[4] Oumayma Bouhamama. "Méthodes numériques pour la résolution du problème inverse en électrocardiographie dans le cas d'anomalies structurelles du tissu cardiaque". PhD thesis. Université de Bordeaux, 2022.

[5] Yves Bourgault, Yves Coudiere, and Charles Pierre. "Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology". In: *Nonlinear analysis: Real world applications* 10.1 (2009), pp. 458-482.

[6] Corify Care. URL: https://corify.es.

[7] Matthijs J M Cluitmans et al. "Noninvasive reconstruction of cardiac electrical activity: update on current methods, applications and challenges". In: *Netherlands Heart Journal* 23 (2015), pp. 301-311.

[8] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. "The heat method for distance computation". In: *Communications of the ACM* 60.11 (2017), pp. 90-99.

[9] Mohamadou Malal Diallo, Yves Coudière, and Rémi Dubois. "A volume source method for solving ECGI inverse problem". In: *International Conference on Functional Imaging and Modeling of the Heart*. Springer. 2021, pp. 551-560.

[10] Josselin Duchateau, Mark Potse, and Remi Dubois. "Spatially coherent activation maps for electrocardiographic imaging". In: *IEEE Transactions on Biomedical Engineering* 64.5 (2016), pp. 1149-1156.

[11] Marc Ethier and Yves Bourgault. "Semi-implicit time-discretization schemes for the bidomain model". In: *SIAM Journal on Numerical Analysis* 46.5 (2008), pp. 2443-2468.

[12] David B Geselowitz. "Description of cardiac sources in anisotropic cardiac muscle: application of bidomain model". In: *Journal of electrocardiology* 25 (1992), pp. 65-67.

[13] David B Geselowitz and WT Miller. "A bidomain model for anisotropic cardiac muscle". In: *Annals of biomedical engineering* 11 (1983), pp. 191-206.

[14] Subham Ghosh and Yoram Rudy. "Application of 11-norm regularization to epicardial potential solution of the inverse electrocardiography problem". In: *Annals of biomedical engineering* 37 (2009), pp. 902-912.

[15] Manuel Haas et al. "Finite element-based space-time total variation-type regularization of the inverse problem in electrocardiographic imaging". In: *CoRR* (2024).

[16] Bin He, Guanglin Li, and Xin Zhang. "Noninvasive imaging of cardiac transmembrane potentials within three-dimensional myocardium by means of a realistic geometry anisotropic heart model". In: *IEEE Transactions on Biomedical Engineering* 50.10 (2003), pp. 1190-1202.

[17] Vitaly Kalinin, Alexander Shlapunov, and Konstantin Ushenin. "On uniqueness theorems for the inverse problem of electrocardiography in the Sobolev spaces". In: *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik* 103.1 (2023), e202100217.

[18] Vitaly Kalinin et al. "On the correctness of the transmembrane potential based inverse problem of ECG". In: *2017 Computing in Cardiology (CinC)*. Ieee. 2017, pp. 1-4.
[19] Amel Karoui. "Méthodes numériques pour la résolution de problèmes inverses en électrocardiographie". PhD thesis. Université de Bordeaux, 2021.

[20] Emma Lagracie. "Reconstruction de la séquence d'activation du coeur par mesure non invasive". PhD thesis. Université de Bordeaux, 2025.

[21] Emma Lagracie et al. "A depth-averaged heart model for the inverse problem of cardiac electrophysiology". In: *Inverse Problems* 41.2 (2025), p. 025002.

[22] GT Lines et al. "Mathematical models and numerical methods for the forward problem in cardiac electrophysiology". In: *Computing and Visualization in Science* 5 (2003), pp. 215-239.

[23] Medtronics. URL: https://www.medtronic.com/en-us/healthcare-professionals/products/cardiac-rhythm/cardiac-mapping/cardioinsight-mapping-system.html.

[24] Colleen C Mitchell and David G Schaeffer. "A two-current model for the dynamics of cardiac membrane". In: *Bulletin of mathematical biology* 65.5 (2003), pp. 767-793.

[25] Mark Potse. "Scalable and accurate ECG simulation for reaction-diffusion models of the human heart". In: *Frontiers in physiology* 9 (2018), p. 370.

[26] Mark Potse et al. "A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart". In: *IEEE Transactions on Biomedical Engineering* 53.12 (2006), pp. 2425-2435.

[27] Andrew J Pullan et al. "The inverse problem of electrocardiography". In: *Comprehensive electrocardiology* 1 (2010), pp. 299-344.

[28] Robert W Roudijk et al. "Comparing non-invasive inverse electrocardiography with invasive endocardial and epicardial electroanatomical mapping during sinus rhythm". In: *Frontiers in physiology* 12 (2021), p. 730736.

[29] Steffen Schuler, Danila Potyagaylo, and Olaf Dössel. "Delay-based regularization for ecg imaging of transmembrane voltages". In: *2019 Computing in Cardiology (CinC)*. Ieee. 2019, Page-1.

[30] Steffen Schuler et al. "Reducing line-of-block artifacts in cardiac activation maps estimated using ECG imaging: A comparison of source models and estimation methods". In: *IEEE Transactions on Biomedical Engineering* 69.6 (2021), pp. 2041-2052.

[31] Guofa Shou, Ling Xia, and Mingfeng Jiang. "Total variation regularization in electrocardiographic mapping". In: *International Conference on Intelligent Computing for Sustainable Energy and Environment*. Springer. 2010, pp. 51-59.

[32] Leslie Tung. "A bi-domain model for describing ischemic myocardial dc potentials." PhD thesis. Massachusetts Institute of Technology, 1978.

[33] Peter M Van Dam et al. "Non-invasive imaging of cardiac activation and recovery". In: *Annals of biomedical engineering* 37 (2009), pp. 1739-1756.

[34] J. Vicente-Puig, J. Chamorro-Servent, and E. et al. Zacur. "Volumetric non-invasive cardiac mapping for accessible global arrhythmia characterization". In: *Commun Med, in press* (2026).

[35] Curtis R Vogel and Mary E Oman. "Iterative methods for total variation denoising". In: *SIAM Journal on Scientific Computing* 17.1 (1996), pp. 227-238.

[36] Jeanne van der Waal et al. "Basis and applicability of noninvasive inverse electrocardiography: a comparison between cardiac source models". In: *Frontiers in physiology* 14 (2023), p. 1295103.

[37] Dafang Wang et al. "Inverse electrocardiographic source localization of ischemia: an optimization framework and finite element solution". In: *Journal of computational physics* 250 (2013), pp. 403-424.

[38] Linwei Wang et al. "Physiological-model-constrained noninvasive reconstruction of volumetric myocardial transmembrane potentials". In: *IEEE Transactions on Biomedical Engineering* 57.2 (2009), pp. 296-315.

[39] Yong Wang and Yoram Rudy. "Application of the method of fundamental solutions to potential-based inverse electrocardiography". In: *Annals of biomedical engineering* 34 (2006), pp. 1272-1288.

[40] Jingjia Xu et al. "Noninvasive transmural electrophysiological imaging based on minimization of total-variation functional". In: *IEEE transactions on medical imaging* 33.9 (2014), pp. 1860-1874.
[41] Yasuo Yamashita and David B Geselowitz. "Source-field relationships for cardiac generators on the heart surface based on their transfer coefficients". In: *IEEE transactions on biomedical engineering* 11 (1985), pp. 964-970.