# Particle Filtering for Non-Deterministic Electrocardiographic Imaging

Emma Lagracie, Luc de Montella

► To cite this version:

Emma Lagracie, Luc de Montella. Particle Filtering for Non-Deterministic Electrocardiographic Imaging. 2025. hal-05272068

HAL Id: hal-05272068

https://inria.hal.science/hal-05272068v1

Preprint submitted on 22 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Particle Filtering for Non-Deterministic Electrocardiographic
Imaging

Emma Lagracie¹,², Luc de Montella²

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

² Inria Bordeaux, F-33400 Talence, France

September 22, 2025

## Abstract

Electrocardiographic imaging (ECGI) aims to non-invasively reconstruct activation maps of the heart from temporal body surface potentials. While most existing approaches rely on inverse and optimization techniques that may yield satisfactory reconstructions, they typically provide a single deterministic solution, overlooking the inherent uncertainty of the problem stemming from its very ill-posed nature, the poor knowledge of biophysical features and the unavoidable presence of noise in the measurements. The Bayesian framework, which naturally incorporates uncertainty while also accounting for temporal correlations across time steps, can be used to address this limitation.

In this work, we propose a low-dimensional representation of the activation sequence that enables the use of particle filtering, a Bayesian filtering method that does not rely on predefined assumptions regarding the shape of the posterior distribution, in contrast to approaches like the Kalman filter. This allows to produce not only activation maps but also probabilistic maps indicating the likelihood of activation at each point on the heart over time, as well as pseudo-probability maps reflecting the likelihood of a point being part of an earliest activation site.

Additionally, we introduce a method to estimate the probability of the presence of a conduction lines of block on the heart surface. Combined with classical reconstruction techniques, this could help discriminate artificial from true lines of block in activation maps. We support our approach with a numerical study based on simulated data, demonstrating the potential of our method.

**Keywords** Inverse problem, Cardiac electrophysiology, ECGi, Bayesian filtering, Particle filtering, Activation probability, Earliest activation site, Line of block

## 1 Introduction

Preceding any contraction, the heart cells undergo an electrical stimulation, called the action potential, which spreads throughout the entire muscle to synchronize the contraction of the heart. If abnormal, this *activation sequence* can trigger malfunctions in the cardiac rhythm. Thus, the electrical activity of the heart provides valuable insight for detecting and identifying cardiac rhythm pathologies.

Electrocardiographic imaging, also called ECGi, aims to non-invasively recover information on the electrical activity of the heart. Using mathematical models of the electrical activity of the heart and torso (source models), ECGi involves solving an inverse problem to reconstruct the bioelectrical sources in the heart (such as transmembrane voltage or extracellular potential) responsible for electrical potential measurements on the torso surface (BSPMs) [1, 2, 3]. In particular, the goal is often to reconstruct the activation sequence of the heart, a quantity of great clinical interest [4].

However, the inverse problem of ECGi is highly ill-posed in the sense of Hadamard, meaning that at least one of the following conditions fails: existence of a solution, uniqueness of the
solution or continuous dependence on the data. In practice, the ill-posedness of the problem is characterized by a very high sensitivity to small changes in the data, making the ECGi problem a very difficult challenge [5]. Moreover, in the clinical setting, source models inevitably involve approximations, such as the omission of heart and torso motion over time, or of the presence of electrical conductivity inhomogeneities caused by organs (lungs, bones, etc.) within the torso. These modeling inaccuracies further increase the complexity of the problem.

Most often, the ECGi problem is formulated as a static Cauchy problem for the Laplace equation in the torso, combined with a Tikhonov regularization term to address its ill-posedness [2, 3, 5, 6]. This allows recovery of the epicardial potential distribution at multiple time instants, from which an activation map can be built [4]. The latter indicates, for each point of the heart, its moment of electrical activation. However, due to the inherent difficulty of the problem, the reconstructions are often not sufficiently accurate. Moreover, in a context involving very partial measurement data, limited knowledge of biophysical electrical properties and strong ill-conditioning, relying on deterministic source models and offering a single deterministic reconstruction does not seem entirely relevant.

Bayesian filtering [7] provides a powerful alternative to tackle the inverse problem of cardiac electrophysiology by rigorously integrating both modeling uncertainties and measurement noise. Furthermore, rather than seeking a single deterministic solution to the ECGi problem, it estimates the full probability distribution of the quantity of interest, called *state* (e.g., extracellular potential or transmembrane voltage), conditioned on the observations (e.g., BSPMs data). This not only avoids forcing a choice among multiple plausible states but also provides tools for confidence assessment. In addition, it naturally introduces state dynamics, which may improve the conditioning of the inverse problem.

Assuming linear Gaussian state dynamics and state-to-observation relationship, the filtering problem can be solved exactly using the Kalman filter. As in related fields for decades [8, 9, 10], the Kalman filter has also been applied to the ECGi problem, see [11, 12, 13] and the comprehensive review [14]. To handle more complex non-linear state dynamics or state-to-observation models, sub-optimal Kalman-based methods, such as the extended Kalman filter and the unscented Kalman filter [15], have been used in electrocardiographic imaging [16, 17]. However, Kalman-based methods still rely on the assumption of a Gaussian posterior state distribution, which can be overly restrictive [16]. In addition, most studies focus only on the maximum a posteriori estimate, thereby identifying the most probable state while discarding valuable information from the full probability density. Although some studies examine how well the filter's covariance predicts the actual reconstruction error [18], restricting the solution to a Gaussian distribution limits error analysis. This also limits the method's ability to deal with multimodal distributions, potentially affecting estimation accuracy [19].

Particle filtering, also known as sequential Monte Carlo methods, provides an alternative to Kalman-based filters without assuming a specific shape for the target distribution [20]. These methods have been successfully applied in many fields [21, 22, 23], and are supported by strong theoretical results [24, 25, 26]. However, they are computationally intensive, which has so far prevented their use in electrocardiographic imaging for reconstructing activation sequences. For instance, the authors in [16] note that particle methods were not feasible for their model due to the high dimensionality of the state (the transmembrane voltage across the heart volume).

To overcome this issue, we introduce a parameterized representation of the cardiac activation sequence, with the goal of drastically reducing the dimensionality of the state. The main parameters used to describe the transmembrane voltage distribution over the whole transmural volume are a set of $l$ points on the heart and $l$ associated positive real numbers, where $l$ is typically less than four. These points can be interpreted as activation centers, or earliest activation sites, and each real number as the radius of a circle (with respect to a geodesic norm adapted to the heart's geometry and conductivity) within which cardiac cells are considered activated. From these parameters and a predefined activation shape, we reconstruct the transmembrane voltage. As a result, the dimensionality of the state (the $l$ centers and $l$ radii) is no longer a barrier to applying particle filtering.
Predefined templates for the transmembrane voltage have already been extensively used in ECGi, as they constitute an efficient way to constrain the space of solutions of the inverse problem [27, 28, 29, 30, 31]. In our setting, it also simplifies the choice of the state evolution model, usually a difficult but essential step in making Bayesian estimation methods effective [11, 32, 33, 16]. In this work, we deliberately choose a very simple one to emphasize that our method can be effective even without a complex dynamic model or strong prior knowledge.

The chosen parameterization makes the state-to-observation relationship highly nonlinear, further motivating the use of particle filtering. More importantly, this particle filtering method enables the computation of clinically relevant quantities beyond the heart's activation map and supports more detailed confidence assessments. For instance, it allows the calculation of the conditional probability, with respect to the data, for a given point on the heart to belong to an earliest activation site or to be activated at a specific time step. To address more complex scenarios, the parameterization can be extended by introducing additional dimensions, such as incorporating a choice among multiple geodesic distance metrics. This extension can be used to assess confidence in the presence of a true conduction block at a specific location. To our knowledge, this is the first approach enabling such an analysis.

The remainder of the paper is organized as follows. In Section 2, we recall the formulation of the electrocardiographic imaging problem, as well as the principles of Bayesian estimation and particle filtering. Section 3 is devoted to describing our low-dimensional representation of the cardiac activation, along with the equations governing its temporal evolution and its relationship to the observed data. The numerical results supporting our method are presented in Section 4, with estimation tools related to the activation sequence detailed in Subsection 4.2, and modifications aimed at distinguishing artificial from true block lines discussed in Subsection 4.3.

# 2 Background

## 2.1 Cardiac Electrophysiology and the ECGI Problem

In mathematical cardiac electrophysiology, the bidomain equations [34, 35] constitute the reference model for computing the dynamics of the electrical potentials and voltages in the heart and torso domains. We denote $\Omega_H$ and $\Omega_T$ the heart and torso domains respectively, $\Omega = \Omega_H \cup \Omega_T$ the whole domain and $\Gamma_T$ the torso surface. The transmembrane voltage $v \in H^1(\Omega_H)$ and the

<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic representation of the domain $\Omega = \Omega_H \cup \Omega_T$.</figcaption></figure>

extracellular and extracardiac potential $u \in H^1(\Omega)$ satisfy the following equations, written in a non-dimensional form: for all $t \in [0, T]$,

$$ 
\left\{
\begin{array}{ll}
\text{div}(\sigma_i \nabla(u+v)) = \partial_t v + f(v,t) & \text{in } \Omega_H, \\
\text{div}((\sigma_i + \sigma_e)\nabla u) = -\text{div}(\sigma_i \nabla v) & \text{in } \Omega_H, \\
\text{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 } \Gamma_T,
\end{array}
\right.
\qquad (1) $$

where $\sigma_{i,e,T}$ are scaled conductivity tensors of the intracellular, extracellular and torso domains respectively, and $f$ is a scaled ionic current, that also depends on a set of coupled ODEs. Due
to the reaction-diffusion equation, during the activation process of the heart, an action potential spreads as a wave across the whole myocardium, offering a very characteristic time course to the transmembrane voltage at each point of the heart. Thus, its shape, contrary to the one of the extracellular potential, can be parametrized with acceptable accuracy using a predefined template.

Based on part of the bidomain equations, and a priori knowledge of cardiac electrophysiology, the inverse ECGi problem consists in recovering the transmembrane voltage $v$ or the extracellular potential $u$ distribution in the heart volume or on the heart surface (epicardium or epi-endocardial surface), to further retrieve an activation map [1, 2, 36, 3]. The most common deterministic approaches seek to reconstruct the epicardial extracellular potential distribution by solving a Cauchy problem for the Laplace equation [2, 37] (static formulation) or use a pre-defined template shape for the transmembrane voltage on the epicardial surface to formulate an activation-time based non-linear inverse problem [27, 37] (intrinsically temporal formulation).

For a clinical use, the quantity of interest is the activation map, in which some features may indicate a pathology. For instance, it is usual to identify the earliest activation sites of the heart, and so the origin of premature ventricular contractions (PVC), or to look at low conduction regions and lines of block in the map, which may trigger ventricular fibrillation [38]. However, we usually lack the means to assess the level of confidence in the reconstructions of the activation maps and its clinical outcome.

## 2.2 Bayesian Estimation

Bayesian estimation is a method for computing the probability distribution of a state, i.e. the quantity of interest such as the extracellular potential or the transmembrane voltage, using observations of that state (e.g. BSPMs data). It consists in sequentially incorporating data through two key models: a state model that describes the stochastic evolution of the state between consecutive measurements, and a measurement equation that characterizes the probabilistic relationship between the state and the observed data.

The state is assumed to be a Markov chain $(X_k)_{k \in \mathbb{N}}$ on $\mathbb{R}^d$, for some $d \in \mathbb{N}^*$. It evolves according to the state equation:

$$X_{k+1} = f_k(X_k, v_k), \qquad (2)$$

where $(v_k)_{k \in \mathbb{N}}$ is a family of independent and identically distributed (i.i.d.) random variables on $\mathbb{R}^d$. These variables model both the random evolution of $X_k$ and uncertainty in the model function $f_k$.

The observation $Y_k \in \mathbb{R}^q$ at time $k$ is linked to the state by the measurement equation:

$$Y_k = h_k(X_k, w_k), \qquad (3)$$

where $(w_k)_{k \in \mathbb{N}}$ is a family of i.i.d. random variables on $\mathbb{R}^q$. These account for both the measurement noise and the uncertainty in the measurement function $h_k$.

Bayesian estimation aims to compute the conditional probability density, denoted $p(X_k | Y_{1:k})$, of the state $X_k$ given the observations $Y_{1:k} := (Y_1, \dots, Y_k)$ up to time $k$. This computation is done iteratively with each new observation, via two steps.

Assuming $p(X_k | Y_{1:k})$ is known, the prediction step forecasts the state density at the next time step using the transition kernel known through the state equation (2):

$$p(X_{k+1} | Y_{1:k}) = \int_{\mathbb{R}^d} p(X_{k+1} | X_k) p(X_k | Y_{1:k}) dX_k. \qquad (4)$$

In the correction step, this prediction is then refined using the new observation and the likelihood known through the measurement equation (3):

$$p(X_{k+1} | Y_{1:k+1}) \propto p(Y_{k+1} | X_{k+1}) p(X_{k+1} | Y_{1:k}). \qquad (5)$$

For example, to apply the Bayesian estimation framework to electrocardiographic imaging, the observation variable $Y_k$ is usually defined as a vector of body surface potentials, and the state
variable $X_k$ can be defined as either the vector of epicardial potentials at each node on the heart surface [39, 11], the transmembrane voltage in the myocardial volume [16], or an updated activation map [17].

Depending on the state one seeks to estimate, choosing a suitable state equation can be more or less complex. The temporal dynamics of the extracellular potential is less explicit than that of the transmembrane voltage, which undergoes an action potential over time. However, this choice of state equation is crucial, as it strongly influences the inverse problem reconstructions. Its deterministic component constrains the temporal evolution of the state via a prior that will have more or less importance depending on the chosen process noise. For example, in [16], Wang and Lui propose to model the dynamics of the transmembrane voltage in the myocardial volume using a monodomain [40, 41] reaction-diffusion model associated with a parameterized ionic model and Gaussian noise. Other works have constrained the extracellular potential on the epicardium to depend on its neighbors at the previous time step [39], in order to better identify the initiation points of the heart's electrical activation. In [17], Lui and He directly estimate the activation map of the myocardium by linking it to torso observations through a prescribed form of the transmembrane voltage and the electrostatic equilibrium part of the bidomain equations (1) linking $v$ to $u$.

The measurement equation is usually simply composed of the state to observation transfer matrix, with an additive Gaussian noise process.

## 2.3 Sequential Importance Resampling

Particle filtering, also known as sequential Monte Carlo, is a set of methods used to approximate the Bayesian estimation equations numerically. It aims to approximate the posterior density $p(X_k | Y_{1:k})$ with a weighted sum of Dirac measures. While many sophisticated particle filtering methods exist and could prove efficient for the problem at hand [42, 43], in this article, we restrict ourselves to the use of the Sequential Importance Resampling (SIR) algorithm [44]. Indeed, our aim is not to compare particle filtering techniques, but rather to propose an efficient way to apply them to ECGi problems. For this purpose, the simplicity of SIR makes it particularly well suited.

The principle of the SIR algorithm is as follows. At time $n = 0$, given a prior distribution density, denoted $p(X_0)$, of the initial state $X_0$, we draw $N$ independent samples $\xi_0^i$, called particles. Each particle is associated with a weight $\omega_0^i$, initially set to $\frac{1}{N}$. The initial density is then approximated by

$$p(X_0) \approx \sum_{i=1}^{N} \omega_0^i \delta_{\xi_0^i}(X_0),$$

where $\delta_{\xi_0^i}$ denotes the Dirac measure centered at $\xi_0^i$.

At the prediction step, the particles evolve according to the state equation (2), and the weights remain unchanged. In practice, we draw $N$ independent realizations $v_k^i$ of the process noise $v_k$, and compute

$$\xi_{k+1}^i = f_k(\xi_k^i, v_k^i), \quad \text{so that} \quad p(X_{k+1} | Y_{1:k}) \approx \sum_{i=1}^{N} \omega_k^i \delta_{\xi_{k+1}^i}(X_{k+1}).$$

At the correction step, the particles remain unchanged, and the weights are updated according to the likelihood of the observation, following equation (5):

$$\omega_{k+1}^i = \frac{\omega_k^i p(Y_{k+1} | \xi_{k+1}^i)}{\sum_{j=1}^{N} \omega_k^j p(Y_{k+1} | \xi_{k+1}^j)}, \quad \text{so that} \quad p(X_{k+1} | Y_{1:k+1}) \approx \sum_{i=1}^{N} \omega_{k+1}^i \delta_{\xi_{k+1}^i}(X_{k+1}).$$

To avoid a situation where only a few particles have non-negligible weights, a phenomenon known as degeneracy, the SIR algorithm includes a resampling step, which is triggered whenever
the effective sample size

$$N_{\text{eff}} = \frac{1}{\sum_{j=1}^{N} (\omega_{k+1}^{j})^{2}}$$

falls below a predefined threshold. A new set of N particles is then drawn, with each particle selected with probability proportional to its weight. This has the effect of duplicating high-weight particles and discarding low-weight ones.

# 3 Method Description

## 3.1 Low-Dimensional State Representation

In this paper, we consider a state variable $X_k$ representing the cardiac activation front. To avoid high-dimensional issues, we leverage the characteristic shape of the activation front [45, 30] and model the electrical activation using $l$ growing balls with respect to a geodesic measure defined on the heart mesh (volume or surface). Each ball is defined by a center and a time-dependent radius, which allows the front to be fully described with only a few parameters. From this activation front, a predefined approximation of the transmembrane voltage is computed at each time step, from which we derive the extracellular and extracardiac potentials. Since the model is based on a low-dimensional state, we can apply particle filtering methods to estimate the centers and radii of the balls, thereby reconstructing the activation map of the heart.

More formally, let us consider, at time $t_k$, $l$ activation centers of the electrical activation of the heart, denoted by $c_k^i \in \mathbb{R}^3$, and referred to as centers. Each activation point is associated with a time-dependent distance $r_k^i \in \mathbb{R}_+$, called the radius. Together, the radii and centers define balls on the heart mesh (volume or surface), with respect to a geodesic distance $d$ defined on that mesh. This distance, which accounts for the shape of the heart, can be weighted according to the conductivity tensors defined on the heart. This, for instance, allows us to consider the fibers of the heart that give a preferential direction of propagation.

The state $X_k$ to be estimated at time $t_k$ is then:

$$X_k = (c_k^1 \cdots c_k^l \quad r_k^1 \cdots r_k^l). \qquad (6)$$

Along with this, we consider a predefined shape of the transmembrane voltage based on a function $V(\xi)$, where $\xi \in \mathbb{R}$ denotes the relative distance between a point on the heart and an activation center. Here, $V$ is expressed as a smoothed Heaviside function:

$$V(\xi) = \begin{cases} 0 & \text{if } \xi < -\text{width}, \\ 1 & \text{if } \xi > +\text{width}, \\ -\frac{\xi^3}{4\text{width}^3} + \frac{3\xi}{4\text{width}} + \frac{1}{2} & \text{if } |\xi| \le \text{width}. \end{cases} \qquad (7)$$

The reconstructed transmembrane voltage at time $t_k$ is then defined as

$$v_{X_k}(x) = \max_{i=1,...,l} V(r_k^i - d(x, c_k^i)). \qquad (8)$$

According to this definition, the transmembrane voltage equals 1 inside each ball defined by a center and a radius (corresponding to the activated region of the heart), and equals 0 outside (corresponding to the resting region). Between these two extreme states, a smooth transition occurs across the activation front. The resting and plateau value of the transmembrane voltage were chosen to be 0 and 1 in accordance to the Mitchell-Shaeffer ionic model [46], but can be adjusted to recover a transmembrane voltage between -90mV and 20mV. It is even possible to consider these extremal values as parameters.

Note that it is not necessary to know the exact number of activation centers but only an upper bound to reconstruct the signal using our method. Indeed, if several balls overlap, the
transmembrane voltage v remains in the activated state 1. Additional centers may also be associated with a null radius, and thus have no impact on the reconstruction.

From the transmembrane voltage v, the extracellular and extracardiac potential u can be computed at each time step using the electrostatic equilibrium equations from the bidomain model (9)

$$ 
\left\{
\begin{array}{ll}
\text{div}((\sigma_i + \sigma_e)\nabla u) = -\text{div}(\sigma_i\nabla v) & \text{in } \Omega_H, \\
\text{div}(\sigma_T \nabla u) = 0 & \text{in } \Omega_T, \\
\sigma_i \nabla (u + v) \cdot n + \sigma_e \nabla u \cdot n = \sigma_T \nabla u \cdot n & \text{on } \partial\Omega_H, \\
\sigma_T \nabla u \cdot n = 0 & \text{on } \Gamma_T,
\end{array}
\right.
\qquad (9) $$

in the case of a volume mesh of the myocardium, or alternatively, the Epicardial Model [47, 48] or the Equivalent Dipole Layer model [49, 50, 51] in the case of a surface epicardial or epi-endocardial mesh.

**Remark** The state representation given in equation (6) will constitute the core of our method and will be used as such in Section 4.2. However, in order to achieve better results, or to more accurately evaluate specific aspects of the activation sequence, it is possible to extend the state. Such an extension is presented in Section 4.3, with the objective of evaluating the probability of the presence of a true conduction block line in the activation map.

## 3.2 State Dynamics and Observation Model

In this paper, we deliberately choose to work with a very simple state equation. Our objective is to test the resilience of our method under minimal modeling assumption. We do not claim however that more elaborate state equations, whether physics-based or data-driven, could not improve the filter's performance, only that we aim to evaluate how well our method performs without relying on such additional inputs.

For the radius, the state equation reads:

$$ r_k^i = \max(r_{k-1}^i + m_{k-1,i}^r, 0). \qquad (10) $$

where $m_{k-1,i}^r \sim \mathcal{N}(0, \Sigma_r^2)$. Note that with this state model, no prior is imposed on the speed of the l activation fronts, and that, for instance, nothing prevent the activation front to go backward. We expect the filter to determine the propagation solely from the data, which allows us to avoid unreliable priors. Note that we choose Gaussian noise model here, but other distributions could also be considered.

Likewise, no prior is imposed on the locations of the centers, which evolve at each time step. Due to the fact that we restrict activation centers to be among heart mesh points, the overall aspect of their temporal evolution take a peculiar form. During the prediction step, we first draw l random variable following an exponential distribution of parameter $\lambda$. Those random variables represent the maximum allowed displacement for each center according to the geodesic metric d defined on the heart mesh. Denoting $\mathcal{E}(\lambda)$, the exponential distribution of parameter $\lambda$, we have:

$$ d_{\max}^i \sim \mathcal{E}(\lambda) \quad \text{and} \quad \mathbf{d}_{\max} = (d_{\max}^1 \cdots d_{\max}^l). $$

The new centers $c_k^i$, for $j = 1, \dots, l$, are then drawn uniformly among the points of the heart mesh lying within the ball of radius $d_{\max,j}$ centered on the previous center $c_{k-1}^i$. Denoting $\mathcal{U}(S)$ the uniform distribution on some set S, the new state equation then reads:

$$ c_k^i \sim \mathcal{U}(B(c_{k-1}^i, d_{\max}^i)) $$

where $B(c_{k-1}^i, d_{\max}^i)$ denotes the ball in $\Omega_H$, with respect to the geodesic distance d, centered on $c_{k-1}^i$ with radius $d_{\max}^i$.

In order to define our observation equation, let O be the linear operator mapping v to u (either the electrostatic equilibrium equations (9) for a volume heart mesh, or the Equivalent
Dipole Layer [50] or Epicardial Model [48] for a surface heart mesh). The observation equation of our filter is written:

$$Y_k = Ov_{X_k}|_{\Gamma_T} + w_k, \quad (11)$$

where $\Gamma_T$ is the torso surface mesh, $|_{\Gamma_T}$ designates the trace operator on this surface and $w_k \sim \mathcal{N}(0, \Sigma_w)$ is a Gaussian measurement noise. Here we use Gaussian observation noise, but other noise models can be considered, provided the probability distribution of $w_k$ can be numerically approximated to compute likelihoods. In our case, the likelihood can be computed as

$$p(Y_k|X_k) \propto \exp \left[ -\frac{1}{2} (Y_k - Ov_{X_k}|_{\Gamma_T})^\top \Sigma_w^{-1} (Y_k - Ov_{X_k}|_{\Gamma_T}) \right]. \quad (12)$$

According to our parametrization of the cardiac activation front, the observation equation is nonlinear, particularly the mapping from $X_k$ to $v_{X_k}$.

Since all electrical potentials are defined up to an additive constant, we also adjust in (12) both the observed data $Y_k$ and the model outputs $Ov_{X_k}|_{\Gamma_T}$ to have zero mean on the torso surface to allow for a meaningful comparison.

Note that when the initial distribution provides little information, the early estimates produced by the particle filter tend to be imprecise. This is because the filter requires sufficient observational data before it can adjust the particle cloud. Since valuable information may be present in these initial moments, this limitation can be a significant drawback. However, nothing prevents us from feeding the observations to the filter in reverse chronological order, and to estimate the evolution of the potential starting from the final time. This approach is expected to allow the filter to gather information more efficiently in the initial stages. To distinguish between these two configurations, we will refer to the standard chronological ordering as the *forward* method (abbreviated as *fwd*), and the reversed ordering as the *backward* method (abbreviated as *bwd*). Without specific mention, the *backward* method is assumed throughout the paper.

The ECGi particle filtering method is summarized in Algorithm 1.

# 4 Numerical Results

## 4.1 Method and Data

### 4.1.1 Data Generation

Data were simulated on a volume heart and torso mesh (17523 myocardial nodes, 2882 torso surface nodes and 44156 total nodes), using the bidomain equations (1) coupled with the Mitchell-Shaeffer ionic model [46]. These equations were discretized spatially using the $P^1$-Lagrange finite elements method, and temporally with a Semi-Implicit Backward Differentiation Scheme of order 2 (SBDF2) [52].

**Propagation test cases** We designed 6 × 2 propagation test cases for the data generation. Six primary test cases were computed with isotropic conductivity tensors $\sigma_i$ and $\sigma_e$ in the heart. Then six twin test cases were calculated with the same parameters but anisotropic conductivity tensors in the heart, using simulated ruled-based fibers from the algorithm of Bayer et al. [53]. In the case of the isotropic (scaled) conductivities, we imposed $\sigma_i = 1I_3$, $\sigma_e = 3I_3$ and $\sigma_T = 2I_3$. In the presence of the ruled-based fibers, we chose the conductivity tensors singular values according to [54]. Denoting by $\ell$ the longitudinal fiber direction and $t$ the transverse direction, we used the following singular values: $\sigma_i^\ell = \sigma_e^\ell = 3$, $\sigma_i^t = 0.3$, $\sigma_e^t = 1.2$ and $\sigma_T = 2I_3$.

Out of the six primary test cases, two contain one initial stimulation site, two contain two initial stimulation sites, with a small delay between the two stimulations, and the last two cases contain one initial stimulation site but also a true conductivity line of block. These conduction line of block were computed by significantly diminishing the intracellular conductivity values inside a designated area. The test cases and their characteristics are summarized in Table 1.
**Algorithm 1 ECGi SIR filter**

**Inputs:**

* The number of particles N and the number l of activation centers per particle.

* The sequence of n body surface potential measurement vectors on the torso: {$Y_k$}$_{k=1}^n$.

**Initialization:**

* Draw $N \times l$ uniformly distributed activation centers {$c_0^{i,j}$ | $i = 1, ..., N$; $j = 1, ..., l$} across the points of the mesh.

* Set the initial radii {$r_0^{i,j}$ | $i = 1, ..., N$; $j = 1, ..., l$} to 1 when using the forward method, or to 150 for the backward method.

* Initialize the particles to $\xi_0^i = (c_0^{i,1}, ..., c_0^{i,l}, r_0^{i,1}, ..., r_0^{i,l})$ for $i = 1, ..., N$.

* Set the weights to $\omega_0^i = \frac{1}{N}$ for $i=1,...,N$.

For $k = 1$ to $n$ do

**Prediction:**

1. Draw new radii and centers for all N particles $\xi_k^i$, $i = 1, ..., N$:

$$r_k^{i,j} = \max(r_{k-1}^{i,j} + \mathcal{N}(0, \Sigma_r^2), 0) \text{ and } c_k^{i,j} \sim \mathcal{U}(B(c_{k-1}^{i,j}, \mathcal{E}(\lambda))).$$

**Correction:**

1. Compute the transmembrane voltage function for each particle:

$$v_{\xi_k^i}: x \mapsto \max_{j=1,...,l} V(r_k^{i,j} - d(x, c_k^{i,j}))$$

2. Compute the unnormalized likelihood of each particle (replace $Y_k$ with $Y_{n-k+1}$ when using the backward method):

$$p(Y_k | \xi_k^i) = \exp \left( -\frac{1}{2} \left( Y_k - Ov_{\xi_k^i|\Gamma_T} \right)^T \Sigma_w^{-1} \left( Y_k - Ov_{\xi_k^i|\Gamma_T} \right) \right)$$

3. Update the weights of the particles: $\omega_k^i = \frac{\omega_{k-1}^i p(Y_k|\xi_k^i)}{\sum_{j=1}^N \omega_{k-1}^j p(Y_k|\xi_k^j)}$

**Resampling:**

If $N_{\text{eff}} = \frac{1}{\sum_{j=1}^N (\omega_k^j)^2} < \frac{N}{3}$ then

1. Resample each particle $\xi_k^i$ independently from the set {$\xi_k^j$}$_{j=1}^N$, with probabilities proportional to the weights {$\omega_k^j$}$_{j=1}^N$.

2. Reset all weights to $\frac{1}{N}$.

end If

end For

**Output:** The posterior density estimate at each time step: $p(X_k | Y_{1:k}) \approx \sum_{i=1}^N \omega_k^i \delta_{\xi_k^i}(X_k)$
<table><thead><tr><th>Stim name</th><th>Line of Block</th><th>Stim localization</th><th>Stim delay</th><th>Isotropic test case name</th><th>Anisotropic test case name</th></tr></thead><tbody><tr><td>Stim 1</td><td>NA</td><td></td><td>NA</td><td>Stim1Iso</td><td>Stim1Ani</td></tr><tr><td>Stim 2</td><td>NA</td><td></td><td>NA</td><td>Stim2Iso</td><td>Stim2Ani</td></tr><tr><td>Delay 1</td><td>NA</td><td></td><td>20 ms</td><td>Delay1Iso</td><td>Delay1Ani</td></tr><tr><td>Delay 2</td><td>NA</td><td></td><td>15 ms</td><td>Delay2Iso</td><td>Delay2Ani</td></tr><tr><td>Block 1</td><td></td><td></td><td>NA</td><td>Block1Iso</td><td>Block1Ani</td></tr><tr><td>Block 2</td><td></td><td></td><td>NA</td><td>Block2Iso</td><td>Block2Ani</td></tr></tbody></table>

<table><caption>Table 1: The 6 × 2 propagation test cases parameters. The first column designates the initial stimulation name. The second column presents, if applicable, the localization of the true conductivity block line. The third column gives a representation of the stimulation on the heart mesh and the fourth column indicates, if applicable, the time delay between the application of different pacings. The two last columns give the associated propagation test case name for the isotropic conductivities as well as for anisotropic conductivities.</caption><thead><tr><th>Stim name</th><th>Line of Block</th><th>Stim localization</th><th>Stim delay</th><th>Isotropic test case name</th><th>Anisotropic test case name</th></tr></thead><tbody><tr><td>Stim 1</td><td>NA</td><td></td><td>NA</td><td>Stim1Iso</td><td>Stim1Ani</td></tr><tr><td>Stim 2</td><td>NA</td><td></td><td>NA</td><td>Stim2Iso</td><td>Stim2Ani</td></tr><tr><td>Delay 1</td><td>NA</td><td></td><td>20 ms</td><td>Delay1Iso</td><td>Delay1Ani</td></tr><tr><td>Delay 2</td><td>NA</td><td></td><td>15 ms</td><td>Delay2Iso</td><td>Delay2Ani</td></tr><tr><td>Block 1</td><td></td><td></td><td>NA</td><td>Block1Iso</td><td>Block1Ani</td></tr><tr><td>Block 2</td><td></td><td></td><td>NA</td><td>Block2Iso</td><td>Block2Ani</td></tr></tbody></table>

**Data interpolation and noise addition** Body surface potential maps (BSPMs) were extracted from the bidomain solutions and interpolated on a second "inverse" torso surface coarser mesh (1186 nodes). Then, the interpolated torso signals were corrupted by adding white noise with a standard deviation equal to 4% of the mean absolute amplitude of the data over time, which roughly corresponds to a signal-to-noise ratio of 30 – 32 dB.

### 4.1.2 Methods for the Filter

**Inverse mesh** In the ECGi procedure, we used a different coarser volume heart and torso mesh of 22114 total nodes, 9174 myocardial nodes and 1186 torso surface nodes. The use of a different mesh than the one used for the data simulation allows to avoid inverse crime.

**Transfer matrix model O** As we used volumetric meshes, the transfer matrix *O* was obtained by discretizing the isotropic electrostatic equilibrium equations (9) with a $P^1$-Lagrange finite elements method. We imposed $\sigma_i = 1I_3$, $\sigma_e = 3I_3$ and $\sigma_T = 2I_3$ for all ECGi resolutions, thus voluntarily introducing a modeling error through *O* when applying the filter on the *Stim1Ani*, *Stim2Ani*, *Delay1Ani*, *Delay2Ani*, *Block1Ani* and *Block2Ani* data.

**Weighted geodesic distance d** The geodesic distances on the heart are computed using a Fast Marching algorithm [55]. To account for tissue conductivities, the weight of each vertex, initially defined as its length, is divided by a correction coefficient. Let *a* and *b* denote the coordinate of the endpoints of the vertex, and let $\sigma_a$ and $\sigma_b$ be the effective conductivity tensors at these points, as defined in the monodomain approximation for proportional intracellular and
extracellular conductivity tensors [40, 34]. The correction coefficient is then given by

$$ \min \left( \sqrt{\frac{a^\top \sigma_a a}{a^\top a}}, \sqrt{\frac{b^\top \sigma_b b}{b^\top b}} \right). \qquad (13) $$

Note that this approach is not a standard way to account for tissue conductivity. However, it was selected for its computational efficiency and to deliberately introduce errors into the information provided to the filtering algorithms, as it would happen in clinical settings. Alternative methods for computing geodesic distance, such as the heat method [56], could also have been employed.

To be consistent with the chosen transfer matrix $O$, we imposed isotropic conductivities $\sigma_i = 1I_3$ and $\sigma_e = 3I_3$ in the geodesic distance calculation. We then expect the distance $d$ to be inaccurate to describe the activation front evolution in the anisotropic test cases *Stim1Ani*, *Stim2Ani*, *Delay1Ani*, *Delay2Ani*, *Block1Ani* and *Block2Ani*. This choice simulates a more realistic situation in which the fiber orientation in the myocardium is generally not known.

**Front width** The front width of the action potential for $V$ was inserted as an a priori value, defined from classical values of the front width in the transmembrane voltage solution of the bidomain model. As shown in previous work [57], some modeling errors in the front width does not generate strong errors on the output body surface potential $u$. Here, the value was set at width = 5mm.

**Covariances and filter parameters** The number of particles is set to $N = 1000$, with each particle containing $l = 3$ centers and radii. Since all test cases contain at most two initial stimulation points, $l$ is always strictly greater than the number of true stimulation points.

The initial locations of the centers are drawn from a uniform distribution over the points of the mesh. For the forward method, the radii $r$ are initialized at 1 mm. With the front width of the potential from (7) set to width = 5 mm, the resulting activated area at initialization is roughly as large as the mesh size on the heart, which corresponds to the smallest detectable value. For the backward method, the radii are set to 150 mm, so that, given the dimensions of the domain, a large part of the heart is activated at the first time step.

Note that with this initialization method, we implicitly assume no prior knowledge of the initial particle distribution. However, if additional information is available, such as a more accurate estimate of the number and location of earliest activation sites, for example from a previous run of the filter, it is possible to adjust the initialization procedure to incorporate this knowledge into the initial particle distribution.

The standard deviation of the white noise governing the evolution of the radii is set to $\Sigma_r = 10$ mm, which is deliberately large to ensure that the filter can follow the signal's evolution. The parameter governing the exponential distribution for the evolution of the centers is set to $\lambda = 5$ mm, to strike a balance between allowing the centers to adjust to the true activation sites and preventing excessive displacement between time steps.

The covariance of the observation equation (11) is set to $\Sigma_w = (0.02)^2 I_d$, where $I_d$ denotes the identity matrix. This value was chosen by estimating the order of magnitude of the Frobenius norm of the difference between the observed data vector and the reconstructed torso potentials from randomly sampled particles. It was then fine-tuned to allow sufficient diversity among particles while avoiding excessive resampling.

## 4.2 Estimation of the Activation Sequence

**Activation maps** We recall that the output of the filter at each time instant $t_k$ is an estimation of the posterior probability density conditioned by the data, given by

$$ p(X_k|Y_{1:k}) \approx \sum_{i=1}^{N} \omega_k^i \delta_{\xi_k^i}(X_k). \qquad (14) $$
Using the notation (8), an estimator of the transmembrane voltage at instant $t_k$ can then be computed as

$$v(x, t_k) \approx \sum_{i=1}^{N} \omega_k^i v_{\xi_k^i}(x). \qquad (15)$$

This estimator enables us to compute activation maps with the maximal spatio-temporal deflec- tion method [4, 58].

<figure><img src="image_3.png" /><figcaption>Figure 2: Reconstructed activation maps for the first four isotropic propagation test cases.</figcaption></figure>

For the test cases *Stim1Iso*, *Stim2Iso*, *Delay1Iso* and *Delay2Iso*, the obtained activation maps are shown in Figure 2 and for their anisotropic counterparts, the recovered activation maps are presented in Figure 3. These visualizations offer an initial glimpse into the potential of the method. When cardiac conductivities are perfectly known in *O* and in the distance *d* (isotropic propagation test cases), the resulting reconstructed maps are of good quality, with a correlation coefficient with the reference greater than 0.97 across all four test cases. In the more realistic test case of inaccurate a priori conductivity tensors in the transfer matrix *O* (i.e. anisotropic propagation test cases), the recovered activation maps are degraded, but key features, such as activation sites, remain clearly identifiable and the correlation coefficient remains above 0.95 for all cases.

Nevertheless, we do not consider the reconstruction of the activation map to be the primary strength of the proposed particle filtering method, therefore we do not include a comparative study with more classical methods on this aspect. Our main objective is to develop nondeter- ministic tools. In the remainder of the paper, we focus on several probabilistic approaches that provide more informative insights than the raw activation maps.

**Activation probability** Using the information provided by the particle filters, it is possible to evaluate, at each time step, the probability that each point of the mesh is activated. This is computed by summing, for each point, the weights of the particles in which the point is considered activated. A point is deemed activated if the transmembrane voltage in that point is superior to a threshold value of 0.5. The probability of being activated at a point *x* of the heart
<figure><img src="image_4.png" /><figcaption>Figure 3: Reconstructed activation maps for the first four anisotropic propagation test cases.</figcaption></figure>

and at time $t_k$ is defined from the $N$ particles of the filter as

$$ \mathbb{P}(v(x, t_k) > 0.5 | Y_{1:k}) \approx \sum_{i=1}^{N} \omega_k^i \mathbb{1}_{[0.5,1]}(v_{\xi_k^i}(x)), \quad (16) $$

where $\mathbb{1}_D$ is the indicator function of the domain $D$ included in $\mathbb{R}$.

Figures 4, 5, and 6 present examples of such activation probability maps for various time steps and test cases. In each case, the activation probabilities computed using either the forward or the backward method are presented side by side with the true transmembrane voltage $v$, ranging between 0 and 1 due to the use of the Mitchell-Shaeffer ionic model.

Several interesting observations can be made. First, when the conductivity tensors in the heart are correctly known by the filter, we obtain fairly accurate results, as shown in Figure 4a. However, as explained earlier, the forward method struggles to produce precise estimates in the early time steps, while the backward method is less accurate toward the end. It can then happen that the probability estimate is incorrect in some regions. For example, still in Figure 4a, the forward method completely misses an activation site at 41 ms. However, in these less accurate phases, the probability maps typically display regions with intermediate activation probabilities. Rather than being a limitation, this reflects the method's capacity to represent and communicate uncertainty. Indeed, it is preferable to have a broad region of uncertainty that contains the true activation zone than to produce a single, deterministic yet misplaced solution.

This uncertainty becomes especially important in test cases where cardiac fibers are present in the true model but omitted in the filter's model. While such a modeling inaccuracy significantly degrades the deterministic activation maps, it has a less severe impact on the activation probability maps, which explicitly express confidence levels for each region. For instance, in Figure 4b, the filter shows high confidence in the initial activation zone, which is accurately estimated, but lower certainty regarding the temporal progression of the activation front.

The study of the Delay test cases provides additional insights. As expected, since the activations occur early in the sequence and are slightly staggered in time, the forward method struggles to correctly identify the activation sites. In contrast, the backward method proves more effective since both activation sites are already active at the final time step, and the filter
<figure><img src="image_5.png" /><figcaption>Figure 4: Probability of being activated at different times for stimulation Stim 2.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Probability of being activated at different times for stimulation *Delay 1*.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Probability of being activated at different times for stimulation *Delay 2*.</figcaption></figure>
can use subsequent observations to better infer the early stages. However, the backward method performs less well in estimating the final stages of the process. This asymmetry underscores the importance of having access to both the forward and backward outputs, which motivates the combined visualization tool presented in the following section.

**Estimation of the earliest activation sites** We now turn our attention to the estimation of the earliest activation sites, rather than the full activation sequence, with the same aim of providing a tool that incorporates confidence information. To the best of our knowledge, no such method has been reported in the literature.

The center estimated by our method, as defined in (6), can be interpreted as a stand-in for these earliest activation sites. To quantify the likelihood of each mesh point to belong to an earliest activation site, we thus compute for each point a pseudo-probability $\tilde{P}(x)$ of being a center by summing the weights of all particles in which the point is selected as a center, for all time steps:

$$ \tilde{P}_k(x|Y_{1:k}) \stackrel{\mathcal{D}}{=} \sum_{i=1}^{N} \omega_k^i \sum_{j=1}^{l} \delta((c_k^j)^i = x), \\ \tilde{P}(x) \stackrel{\mathcal{D}}{=} \sum_{k=1}^{n} \tilde{P}_k(x|Y_{1:k}) \qquad (17) $$

where $(c_k^j)^i$ denotes the j-th center of the i-th particle at time $t_k$.

It is important to emphasize that this quantity is not a true probability, mostly because the values can exceed 1, as they range from 0 to the maximum number of centers represented in the state times the number of time steps. This is due to the fact that the actual number of distinct earliest activation sites is unknown, making it difficult to normalize these values into proper probabilities. Indeed multiple centers within a single particle may correspond to the same physical earliest activation site, and the number of time a center is present in a particle has no influence on the particle's likelihood.

Figures 7 and 8 show the maps of these pseudo-probabilities, obtained by averaging the results of the forward (index fwd) and backward (index bwd) methods, for the first four isotropic test cases and the first four anisotropic test cases respectively. On each map, the true stimulation sites are indicated by black dots.

When the filter uses the exact electrical conductivity tensors (isotropic test cases, Figure 7), the earliest activation sites are located with high accuracy. In the test cases *Stim1Iso*, *Stim2Iso* and *Delay2Iso*, the true activation sites are even identified in a unique manner, meaning that no other possible stimulation sites appear on the map. In the *Delay1Iso* test case, where one of the two true stimulation sites is located near the right ventricular outflow tract, a hardly observable area of the heart, the two correct stimulation sites are identified among four possible sites. Near the right ventricular outflow tract, the potential stimulation zone appears larger, with a lower, more spread pseudo-probability. These results are encouraging as we expect the filter to display all the likely possibilities, and not only one deterministic guess. Moreover, note that the temporal probability of being activated shown in Figure 5a can also allow to discriminate true from false activation sites with a higher level of confidence.

In contrast, when significant modeling errors are introduced through the absence of fibers in the filter heart model (anisotropic data test cases, Figure 8), the estimates are degraded, and the predicted centers tend to be slightly misplaced. However, this error is balanced by a lower confidence in the prediction and a correspondingly broader uncertainty region. Except in the cases where a stimulation occurs near the right ventricular outflow tract, where signal reconstruction proves particularly difficult (*Stim2Ani* and *Delay1Ani*), the true activation sites remain contained within the uncertainty region. This demonstrates the filter's ability to effectively represent and communicate uncertainty in challenging settings.

In summary, the combined multiple visualization tools our filter offers may allow to refine our vision and understanding of the reconstructed output of ECGi.
<figure><img src="image_8.png" /><figcaption>Figure 7: Pseudo-probability of belonging to an initial activation site for the first four isotropic propagation test cases.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Pseudo-probability of belonging to an initial activation site for the first four anisotropic propagation test cases.</figcaption></figure>
4.3 Distinguishing artificial from real conduction lines of block

<figure><img src="image_10.png" /><figcaption>Figure 9: Activation maps obtained with only one homogeneous geodesic distance. The conduc-
tion lines of block cannot be reconstructed.</figcaption></figure>

If we now look at the last 2 × 2 propagation test cases *Block 1* and *Block 2*, which include a real conduction line of block (see Table 1), we expect the filter not to be able to recover the lines of block in the activation maps. Indeed, our chosen parametrization of the activation front shape, trough balls in the isotropic and homogeneous geodesic distance, prevents from recovering a non circular activation shape in the myocardium. This is clearly observed in the recovered activation maps presented in Figure 9 for the *Block1Iso*, *Block1Ani*, *Block2Iso* and *Block2Ani* test cases.

Nevertheless, recovering true regions of low conduction velocity is of interest to clinicians as it
can indicate pathology. Furthermore, it is well known that the classical Cauchy-Tikhonov-based
approache to ECGi tends to produce too many lines of block in the reconstructed activation
maps, among which some are true, and many are artifacts [4, 38]. It would thus be interesting
to overcome the intrinsic limitation of our particle filtering model, while offering a tool for
discriminating true from artificial conduction lines of block in the activation maps.

We thus propose to adapt our method, not to detect lines of block directly, but to help
distinguish between artificial and true lines of block. Assuming that a pre-existing activation
map has been computed using any method (e.g. solving the Cauchy problem for the Laplace
equation), and that lines of block have been identified on this map, it is possible to compute
new geodesic distances and operator *O* that account for the underlying changes in conductivity.
In that case, electrical conductivities in (13) and in *O* can be tuned to produce a distance map
containing a line of block similar, though not perfectly identical, to the one in the activation
map.

We then adapt the filtering method by extending the state defined in (6) with a discrete
component. This discrete state allows each particle to choose from a predefined set of geodesic
distances and corresponding transfer matrix *O*. The distribution of particles across the different
geodesic metrics, those with or without a line of block for instance, provides an evaluation of
the confidence in the actual presence of such a conduction block in the heart. Denoting by $m_k^i$
the mode, composed of a geodesic distance and its associated transfer matrix, of particle *i* at
time $t_k$, the probability reconstructed by the filter that a given geodesic distance $d$ is valid at time $t_k$ is given by:

$$p_k(d) = \sum_{i=1}^{N} \omega_k^i \delta_{m_k^i}(d). \qquad (18)$$

In practice, each particle is initialized by uniformly selecting one of the available geodesic distances. Then, at each time step, each particle has a 99% chance of keeping the same distance. Otherwise, it uniformly selects a new distance from the available options. Allowing particles to switch modes not only enables them to correct themselves over time, but also allows the filter to test the presence of multiple candidate lines of block, even if no single geodesic metric accounts for all of them simultaneously. Indeed, the filter's estimate can switch between distances based on which one best fits the data at each time step. In addition to this mode evolution, the resampling steps tend to duplicate particles associated with the most probable distance metric and corresponding operator $O$.

<figure><img src="image_11.png" /><figcaption>Figure 10: Example of the computed geodesic distance from the points in the red zone, with an a priori block. The reconstruction method of equation (13), intentionally imperfect, produces the block's atypical shape.</figcaption></figure>

Using approximations of the geodesic distances based on the *Block 1* and *Block 2* test cases, we ran 50 *forward* and 50 *backward* filters on the 2×2 test cases that include a line of block. Since each run is independent, the computation can be fully parallelized, avoiding any additional computation time. Each filter was provided with the two distance maps derived from the *Block 1* and *Block 2* test cases, as well as the isotropic distance. We thus expect the filter to be able to have a high confidence in the presence of the true line of block, but also to discard the false one. Figure 11 shows, for each test case, the mode distribution, computed using (18), averaged across all filters. As for earliest activation estimate, computing an average across the *forward* and *backward* method allows us to combine the strengths of those methods.

For test cases based on *Block 2*, regardless of the modeling error due to the presence of fibers, the filters show a very high confidence in the presence of the real line of block, as illustrated in Figures 11a and 11b. In contrast, the block from *Block 1* test cases is more difficult to detect with high confidence. For filters that do not make an assumption error regarding the presence of fibers, the filter is confident in the presence of the real block during the first 60 ms, as shown in Figure 11a. Beyond this point, all geodesic distances exhibit roughly the same likelihood. According to the activation maps in Figure 9, the 60 ms mark corresponds to the moment when the block is crossed, after which its effect becomes very weak. As the confidence level in the first part is lower than that of the filter associated with the *Block 2* test cases, this explains the indecision of the filter in the latter part of the simulation. For test case *Block1Ani*, due to
<figure><img src="image_12.png" /><figcaption>Figure 11: Evolution of the probability distribution over geodesic distance modes, averaged across 50 *forward* and 50 *backward* filter outputs, for test cases containing a true block line. The blue line indicates the actual block line present in the heart for each case.</figcaption></figure>

the lower observability of the block combined with the modeling error caused by the omission of fibers, the filters tend to overlook the real block and instead show confidence in the absence of a block, as illustrated in Figure 11c. Note, however, that throughout the various test cases, the filters have consistently been able to discard the false block, regardless of the modeling errors.

Activation maps can also be recomputed, using the output of the extended-state filter. Such maps are shown in Figure 12. With the extended state and choice between three possible distance maps and operators *O*, we are able to reconstruct accurate activation maps, containing the true conduction lines of block for all four test cases *Block1Iso*, *Block1Ani*, *Block2Iso* and *Block2Ani*. For the *Block1Ani* test case, which mode distribution across time did not indicate the presence of a true line of block (see Figure 11c), we are still able to observe a reconstructed line of block in the activation map, though smaller than in the reference map. For the *Block2Iso* and *Block2Ani* cases, the line of block localization and presence is well recovered. Some singularities appear on both sides of the discontinuity line. These are in fact caused by approximations in the reconstructed geodesic distance, which, for instance, may exaggerate the strength of the block line compared to the real one. Figure 10 provides a glimpse of these imperfections in the geodesic distance computed using formula (13).

These results tend to validate the ability of the extended filter to discriminate true from artificial lines of block, even with imprecision in the estimated distance matrix and conductivity tensors in *O*.

# 5 Conclusion

In this paper, we introduced a new method for estimating various characteristics of the cardiac activation sequence from body surface potential maps. This approach is based on particle filtering techniques, which offer several advantages. First, it allows for more complex transition and observation models than the standard linear Gaussian assumptions, making it possible to
<figure><img src="image_13.png" /><figcaption>Figure 12: Activation maps obtained with an extended state which includes the choice between three geodesic distances: one homogeneous, one heterogeneous containing the true a priori conduction line of block and one heterogeneous containing a false a priori conduction line of block. The true conduction lines of block are recovered in the maps or the four test cases.</figcaption></figure>

incorporate more realistic noise distributions. Second, it avoids the need to assume a predefined shape for the posterior distribution, as is required by Kalman-based filters. This added flexibility enables a more nuanced evaluation of the confidence in the estimates, beyond a simple covariance matrix.

In addition to the classical activation map, we proposed three new estimators based on probabilistic interpretations. The first provides, for each point and time step, the probability of being activated at that time. This estimator offers a detailed view of the activation process, enabling the identification of activation sites and activation times, along with an assessment of the confidence in these estimates. To offer a simpler tool for identifying the earliest activation sites and evaluating their reliability, we introduced a map that assigns to each point a pseudo-probability of belonging to an earliest activation site. Finally, to complement existing deterministic approaches, we proposed a method for assessing whether a detected line of block is likely to be physiological or an artifact, thereby introducing a probabilistic dimension to otherwise deterministic frameworks.

We tested our method on several simulated test cases using a myocardial volume mesh, though the method could equally be applied to a surface mesh of the heart. In the most optimistic cases, where the conductivities of the heart tissue were known, the method was able to recover the correct number of earliest activation sites, through the temporal activation probability maps (Figures 4, 5, 6). At the same time, the true earliest activation sites lay within high pseudo-probability regions of the map indicating the likelihood of belonging to an earliest activation site (Figure 7). In fact, rather than being opposed, these different views can be combined to improve the interpretation of the results, for example by discarding incoherent earliest activation sites. To better mimic real-world situations, where these conductivities are typically unknown, we also considered cases with fibers present in the heart during data generation but not accounted for in the filter. Even with the introduction of this model mismatch, our method still produced satisfactory results. Although the estimates were often shifted, they generally remained within the computed uncertainty regions, whether for the activation time at a given point or for the
earliest activation sites. Moreover, the method proved effective in estimating the number of
earliest activation sites.

Due to the imposed shape of the activation front, made of overlapping activation balls, our method still faces some limitations. For example, if a line of block is present in the heart but not modeled in the filter, the method fails to detect it. To address this issue, we proposed an extension that successfully overcomes the problem. Other challenging issues will likely arise in practical applications, but similar state extensions could be employed to address them effectively.

Though we voluntarily introduced modeling errors in the filter model, we still used simulated data, which are well known to produce ECGi reconstructions of higher quality than clinical or experimental data. Indeed, in a clinical setting, a lot of unknowns, modeling errors and noise, add to the difficulty of the ECGi problem. Moreover, real word data are usually obtained through 148 to 255 electrodes vests, which is way less than the ~ 1000 body surface data points of our torso surface mesh. Our methodology clearly extends to the usual torso electrode vest, provided a coarse mesh of the torso is available, but needs to be assessed in such a setup. Thus future work should focus on evaluating our new method on wider and more challenging experimental datasets.

It should be noted that certain a priori information can also be added to the filter in order
to improve its performance. One could imagine introducing a priori fiber directions into the
transfer matrix model *O* and into the geodesic distances, or using a prior deterministic ECGi
reconstruction to deduce the initialisation state *X*₀.

Acknowledgments

The authors would like to thank Lisl Weynans and Yves Coudière for taking the time to proofread
the paper.

References

[1] Andrew J Pullan, Leo K Cheng, Martyn P Nash, Alireza Ghodrati, Rob MacLeod, Dana H Brooks, et al. The inverse problem of electrocardiography. *Comprehensive electrocardiology*, 1:299–344, 2010.

[2] Laura Bear, Phillip S Cuculich, Olivier Bernus, Igor Efimov, and Rémi Dubois. Introduction to noninvasive cardiac mapping. *Cardiac electrophysiology clinics*, 7(1):1–16, 2015.

[3] Matthijs J M Cluitmans, RLM Peeters, RL Westra, and PGA Volders. Noninvasive reconstruction of cardiac electrical activity: update on current methods, applications and challenges. *Netherlands Heart Journal*, 23:301–311, 2015.

[4] Steffen Schuler, Matthias Schaufelberger, Laura R Bear, Jake A Bergquist, Matthijs JM Cluitmans, Jaume Coll-Font, Önder N Onak, Brian Zenger, Axel Loewe, Rob S MacLeod, et al. Reducing line-of-block artifacts in cardiac activation maps estimated using ecg imaging: A comparison of source models and estimation methods. *IEEE Transactions on Biomedical Engineering*, 69(6):2041–2052, 2021.

[5] Faker Ben Belgacem. Why is the cauchy problem severely ill-posed? *Inverse problems*, 23(2):823, 2007.

[6] AN Tikhonov and V Arsenin. Solutions of ill-posed problems. (No Title), 1977.

[7] Simo Särkkä and Lennart Svensson. *Bayesian filtering and smoothing*, volume 17. Cambridge university press, 2023.

[8] J.W. Phillips, R.M. Leahy, and J.C. Mosher. Meg-based imaging of focal neuronal current sources. *IEEE Transactions on Medical Imaging*, 16(3):338–348, 1997.
[9] S. Baillet and L. Garnero. A bayesian approach to introducing anatomo-functional priors in the eeg/meg inverse problem. *IEEE Transactions on Biomedical Engineering*, 44(5):374-385, 1997.

[10] David M Schmidt, John S George, and Chris C Wood. Bayesian inference applied to the electromagnetic inverse problem. *Human brain mapping*, 7(3):195-212, 1999.

[11] Umit Aydin and Yesim Serinagaoglu Dogrusoz. A kalman filter-based approach to reduce the effects of geometric errors and the measurement noise in the inverse ecg problem. *Medical & Biological Engineering & Computing*, 49(9):1003-1013, 2011.

[12] Yves Goussard, Diane Joly, and Pierre Savard. Time-recursive solution to the inverse problem of electrocardiography. *IEEE Transactions on Biomedical Engineering*, 1993.

[13] Walther HW Schulze, Francesc Elies Henar, Danila Potyagaylo, Axel Loewe, Matti Stenroos, and Olaf Dössel. Kalman filter with augmented measurement model: an ecg imaging simulation study. In *International Conference on Functional Imaging and Modeling of the Heart*, pages 200-207. Springer, 2013.

[14] Yesim Serinagaoglu Dogrusoz. Statistical estimation applied to electrocardiographic imaging. In *2019 12th International Conference on Measurement*, pages 2-9. IEEE, 2019.

[15] Eric A Wan and Rudolph Van Der Merwe. The unscented kalman filter for nonlinear estimation. In *Proceedings of the IEEE 2000 adaptive systems for signal processing, communications, and control symposium (Cat. No. 00EX373)*, pages 153-158. Ieee, 2000.

[16] Linwei Wang, Heye Zhang, Ken CL Wong, Huafeng Liu, and Pengcheng Shi. Physiological-model-constrained noninvasive reconstruction of volumetric myocardial transmembrane potentials. *IEEE Transactions on Biomedical Engineering*, 57(2):296-315, 2009.

[17] Chenguang Liu and Bin He. Noninvasive estimation of global activation sequence using the extended kalman filter. *IEEE Transactions on Biomedical Engineering*, 58(3):541-549, 2010.

[18] Yesim Serinagaoglu, Dana H Brooks, and Robert S MacLeod. Bayesian solutions and performance analysis in bioelectric inverse problems. *IEEE Transactions on Biomedical Engineering*, 52(6):1009-1020, 2005.

[19] Namrata Vaswani. Particle filtering for large-dimensional state spaces with multimodal observation likelihoods. *IEEE Transactions on Signal Processing*, 56(10):4583-4597, 2008.

[20] Jun S Liu and Rong Chen. Sequential monte carlo methods for dynamic systems. *Journal of the American statistical association*, 93(443):1032-1044, 1998.

[21] Randy L Haupt and Douglas H Werner. Genetic algorithms in electromagnetics. John Wiley & Sons, 2007.

[22] Hedibert F Lopes and Ruey S Tsay. Particle filters and bayesian inference in financial econometrics. *Journal of Forecasting*, 30(1):168-209, 2011.

[23] Enzo Iglésis and Luc de Montella. Target motion analysis using angular measurements and underwater sound pressure levels. working paper or preprint, 2025.

[24] Pierre Del Moral and Alice Guionnet. On the stability of measure valued processes with applications to filtering. *Comptes Rendus de l'Académie des Sciences - Series I - Mathematics*, 329(5):429-434, 1999.

[25] Nicolas Chopin, Omiros Papaspiliopoulos, et al. *An introduction to sequential Monte Carlo*, volume 4. Springer, 2020.
[26] Michel Caffarel, Pierre Del Moral, and Luc de Montella. On the mathematical foundations of diffusion monte carlo. *Journal of Mathematical Physics*, 66(1), 2025.

[27] Peter M Van Dam, Thom F Oostendorp, André C Linnenbank, and Adriaan Van Oosterom. Non-invasive imaging of cardiac activation and recovery. *Annals of biomedical engineering*, 37:1739–1756, 2009.

[28] Gwladys Ravon, Rémi Dubois, Yves Coudière, and Mark Potse. A parameter optimization to solve the inverse problem in electrocardiography. In *International Conference on Functional Imaging and Modeling of the Heart*, pages 219–229. Springer, 2017.

[29] Jérôme Fehrenbach and Lisl Weynans. Source and metric estimation in the eikonal equation using optimization on a manifold. *Inverse Problems and Imaging*, 17(2):419–440, 2023.

[30] Piero Colli Franzone, Luciano Guerri, and Stefania Tentoni. Mathematical modeling of the excitation process in myocardial tissue: Influence of fiber rotation on wavefront propagation and potential field. *Mathematical biosciences*, 101(2):155–235, 1990.

[31] Thomas Grandits, Alexander Effland, Thomas Pock, Rolf Krause, Gernot Plank, and Simone Pezzuto. Geasi: Geodesic-based earliest activation sites identification in cardiac models. *International journal for numerical methods in biomedical engineering*, 37(8):e3505, 2021.

[32] Taha Erenler and Yesim Serinagaoglu Dogrusoz. Ml and map estimation of parameters for the kalman filter and smoother applied to electrocardiographic imaging. *Medical & Biological Engineering & Computing*, 57:2093–2113, 2019.

[33] Yesim Serinagaoglu Dogrusoz and Taha Erenler. Use of simulated data for the estimation of prior models in kalman filter-based ecgi. In *2020 Computing in Cardiology*, pages 1–4. IEEE, 2020.

[34] Glenn Terje Lines, Per Grottum, and Aslak Tveito. Modeling the electrical activity of the heart: a bidomain model of the ventricles embedded in a torso. *Computing and Visualization in Science*, 5:195–213, 2003.

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

[36] Lei Li, Julia Camps, Blanca Rodriguez, and Vicente Grau. Solving the inverse problem of electrocardiography for cardiac digital twins: A survey. *IEEE Reviews in Biomedical Engineering*, 2024.

[37] Jeanne van der Waal, Veronique Meijborg, Ruben Coronel, Rémi Dubois, and Thom Oostendorp. Basis and applicability of noninvasive inverse electrocardiography: a comparison between cardiac source models. *Frontiers in physiology*, 14:1295103, 2023.

[38] Josselin Duchateau, Frédéric Sacher, Thomas Pambrun, Nicolas Derval, Judit Chamorro-Servent, Arnaud Denis, Sylvain Ploux, Mélèze Hocini, Pierre Jaïs, Olivier Bernus, et al. Performance and limitations of noninvasive cardiac activation mapping. *Heart rhythm*, 16(3):435–442, 2019.

[39] Shijie Zhou, John L Sapp, Fady Dawoud, and B Milan Horáček. Localization of activation origin on patient-specific epicardial surface by empirical bayesian method. *IEEE Transactions on Biomedical Engineering*, 66(5):1380–1389, 2018.

[40] Mark Potse, Bruno Dubé, Jacques Richer, Alain Vinet, and Ramesh M Gulrajani. A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart. *IEEE Transactions on Biomedical Engineering*, 53(12):2425–2435, 2006.
[41] GT Lines, ML Buist, P Grottum, AJ Pullan, J Sundnes, and A Tveito. Mathematical models and numerical methods for the forward problem in cardiac electrophysiology. *Computing and Visualization in Science*, 5:215-239, 2003.

[42] Branko Ristic, Sanjeev Arulampalam, and Neil Gordon. *Beyond the Kalman filter: Particle filters for tracking applications*. Artech house, 2003.

[43] Jos Elfring, Elena Torta, and René Van De Molengraft. Particle filters: A hands-on tutorial. *Sensors*, 21(2):438, 2021.

[44] Neil J Gordon, David J Salmond, and Adrian FM Smith. Novel approach to nonlinear/non-gaussian bayesian state estimation. In *IEE proceedings F (radar and signal processing)*, volume 140, pages 107-113. IET, 1993.

[45] James P Keener. An eikonal-curvature equation for action potential propagation in myocardium. *Journal of mathematical biology*, 29(7):629-651, 1991.

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

[47] Emma Laggracie, Yves Bourgault, Yves Coudière, and Lisl Weynans. A depth-averaged heart model for the inverse problem of cardiac electrophysiology. *Inverse Problems*, 41(2):025002, 2025.

[48] Emma Laggracie, Lisl Weynans, and Yves Coudière. Effects of the insertion of epicardial anisotropy versus isotropy in the inverse problem of electrocardiography. In *FIMH 2025-13th Functional Imaging and Modeling of the Heart International Conference*, number 15672, pages 83-93. Springer Nature Switzerland, 2025.

[49] Yasuo Yamashita and David B Geselowitz. Source-field relationships for cardiac generators on the heart surface based on their transfer coefficients. *IEEE transactions on biomedical engineering*, (11):964-970, 1985.

[50] David B Geselowitz. Description of cardiac sources in anisotropic cardiac muscle: application of bidomain model. *Journal of electrocardiology*, 25:65-67, 1992.

[51] Jeanne van der Waal, Veronique Meijborg, Ruben Coronel, Rémi Dubois, and Thom Oostendorp. Basis and applicability of noninvasive inverse electrocardiography: a comparison between cardiac source models. *Frontiers in Physiology*, Volume 14 - 2023, 2023.

[52] Marc Ethier and Yves Bourgault. Semi-implicit time-discretization schemes for the bidomain model. *SIAM Journal on Numerical Analysis*, 46(5):2443-2468, 2008.

[53] Jason D Bayer, Robert C Blake, Gernot Plank, and Natalia A Trayanova. A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. *Annals of biomedical engineering*, 40:2243-2254, 2012.

[54] Mark Potse. Scalable and accurate ecg simulation for reaction-diffusion models of the human heart. *Frontiers in physiology*, 9:370, 2018.

[55] J A Sethian. A fast marching level set method for monotonically advancing fronts. *Proceedings of the National Academy of Sciences*, 93(4):1591-1595, 1996.

[56] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. The heat method for distance computation. *Communications of the ACM*, 60(11):90-99, 2017.

[57] Emma Laggracie, Lisl Weynans, and Yves Coudière. Comparison of two formulations for computing body surface potential maps. In *2023 Computing in Cardiology (CinC)*, volume 50, pages 1-4. IEEE, 2023.
[58] Emma Lagracie, Lisl Weynans, and Yves Coudière. Assessment of a threshold method for computing activation maps from reconstructed transmembrane voltages. In *51st international Computing in Cardiology conference-CinC 2024*, 2024.