GEASI: Geodesic-based Earliest Activation Sites Identification in cardiac models

The identification of the initial ventricular activation sequence is a critical step for the correct personalization of patient-specific cardiac models. In healthy conditions, the Purkinje network is the main source of the electrical activation, but under pathological conditions the so-called earliest activation sites (EASs) are possibly sparser and more localized. Yet, their number, location and timing may not be easily inferred from remote recordings, such as the epicardial activation or the 12-lead electrocardiogram (ECG), due to the underlying complexity of the model. In this work, we introduce GEASI (Geodesic-based Earliest Activation Sites Identification) as a novel approach to simultaneously identify all EASs. To this end, we start from the anisotropic eikonal equation modeling cardiac electrical activation and exploit its Hamilton--Jacobi formulation to minimize a given objective function, e.g. the quadratic mismatch to given activation measurements. This versatile approach can be extended to estimate the number of activation sites by means of the topological gradient, or fitting a given ECG. We conducted various experiments in 2D and 3D for in-silico models and an in-vivo intracardiac recording collected from a patient undergoing cardiac resynchronization therapy. The results demonstrate the clinical applicability of GEASI for potential future personalized models and clinical intervention.


Introduction
In this work, we address the central question of identifying earliest activation sites (EASs) in a propagation model for ventricular activation. In a healthy (human) heart, ventricles are activated via a specific pathway that originates in the atrio-ventricular node, continues in the His bundle and the Purkinje network, to eventually spread in the myocardium through Purkinje-myocardial junctions [1]. These junction points can be effectively modeled by a discrete set of EASs, that form the initial condition of the propagation model. Unfortunately, the precise structure of the set of EASs (defined by their number, location and timing) cannot be detected in vivo, and rule-based approaches are limited by inter-patient variability. More importantly, severe pathological conditions such as intraventricular conduction disorders are directly associated with partially malfunctioning activation pathways, hence corresponding to a pathological set of EASs. A correct and possibly automatic identification of EASs from non-invasive or minimally-invasive recordings is therefore of high clinical relevance, especially in selecting the optimal treatment for the patient [2].
A particularly suitable propagation model in the context of EASs is the anisotropic eikonal equation, which was originally exploited as a convenient approximation of the monodomain and bidomain models [3,4], but is nowadays more often utilized for its computational efficiency [5,6]. This work, however, leverages the eikonal model from a novel perspective, based on Hamilton-Jacobi formalism and geodesics [7], to enable a gradient-based approach for localizing the EASs termed GEASI (Geodesic-based Earliest Activation Sites Identification). In detail, we start from the anisotropic eikonal equation as a common model for cardiac electrophysiology [8], in which the EASs define boundary conditions at specific sites. For numerical reasons, the eikonal equation is solved using the Fast Iterative Method (FIM) [9]. The main goal of our approach is the minimization of a given objective functional depending on the solution of the anisotropic eikonal equation as a function of the EASs. Here, a feasible optimization strategy involves the Hamilton-Jacobi formalism, which promotes a tractable derivative with respect to the EASs [7]. Note that this derivative is geometrically related to the tangent of the geodesic at the EASs. In this respect, a geodesic connects an EAS such as a Purkinje entry point to an observation through a path of minimum distance in a predefined metric. Finally, we exploit the aforementioned methods to introduce GEASI, which in its core employs a quadratic mismatch between the eikonal solution and the measurements, such as observations of an electro-anatomical mapping, in the objective function. In summary, GEASI can therefore fit the parameters of the eikonal model to clinical recordings in a very efficient and flexible manner.
We emphasize that GEASI is not limited to this quadratic objective function and can straightforwardly be extended to other scenarios. In this work, we additionally investigate two such extensions: the topological gradient designed for estimating the number of EASs and fitting of EASs of an eikonal-based ECG model to a clinically recorded ECG.
Changing the number of EASs and the effect of this action on the objective function is evaluated by means of the topological gradient. The concept of topological gradients can be readily introduced via the Hamilton-Jacobi theory. Here, we consider the splitting of a single EAS into a pair of two EASs symmetrically arranged at infinitesimal distance along a given direction in a GEASI is able to identify these parameters by exploiting the Hamilton-Jacobi formulation. This formulation allows for an efficient optimization scheme for minimizing mismatches to either activation maps or resulting quantities such as the ECG. Using the topological gradient we can additionally estimate the number of EASs.
dipole-like fashion. Thus, the topological gradient provides a criterion to decide whether an EAS should be split or not. In particular, this approach promotes a simple model as a starting point with too little complexity to represent the measurement data and increase the number of source sites until the encountered activations are properly approximated.
In combination with the pseudo-bidomain model, a template-based action potential and the lead field theory, the eikonal model also results in an almostreal-time ECG simulator [5,6] with remarkable physiological accuracy [8]. Here, we extend a previously proposed approach [10], based on this ECG model, to solve the inverse ECG problem, i.e. we exploit GEASI to localize EASs purely from ECG data. Numerical experiments in Section 5 demonstrate that the proposed approach is capable of finding the optimal EASs even in high-fidelity cardiac models. We visually summarized GEASI and its applications in Figure 1.

Related Work
In what follows, we briefly review similar and related approaches to GEASI.
From a physiological perspective, EASs can be derived from an automatically generated Purkinje network [11], closely following the actual anatomy of the heart. The approach, anatomically-tailored but not patient-specific, is indicated in the case of a generic healthy activation and a (complete or partial) bundle branch block. When dense endocardial mapping data are available, the Purkinje network can be estimated automatically [12,13,14]. In the method proposed by Palamara et al., the Purkinje network is created from intra-cardiac measurements by dividing the endocardium into regions of influence for each Purkinje entry point. For this purpose, an isotropic eikonal equation is solved for each measurement point to compute the regions of influence in a Voronoidiagram like fashion. According to fractal laws, Purkinje entry points are subsequently either moved, deleted or generated to better fit the observed activation on the endocardium. The connection with GEASI becomes apparent once we consider the underlying problem in terms of geodesics and regions of influence, further discussed in Section 2. We can observe that GEASI is a generalization of the aforementioned approach, since it gives rise to Voronoi partitions with nonlinear boundaries. Thus, GEASI can be applied to heterogeneous conduction velocities and heterogeneous fiber directions.
An alternative formulation of the Purkinje network is based on a very sparse set of EASs embedded in thin, fast-conducting layers in both ventricular endocardia [15]. In some sense, this approach can be referred to as lumped Purkinje network formulation. Here, the number of EASs is drastically reduced. In fact, a few sites per chamber are generally sufficient to correctly capture the activation and reproduce the surface ECG [16]. Overall, the parameters of the eikonal model with lumped Purkinje network are just the location, the number and the activation onset of the EASs, as well as the conduction velocities in the thin layers and the myocardium. The problem of fitting these parameters to clinical data has already been considered in the literature for the conduction velocities [17,18,19,20]. However, the optimization of EASs received limited attention so far with only a few works dealing with activation onsets [21] or locations [22]. The simultaneous optimization of EASs (especially their number) and conduction velocity has been analyzed only very recently in [23,20,10].
In the work by Kunisch et al., the authors recast the problem of localizing EASs as a shape optimization problem. In their viscous eikonal formulation, the EASs are modelled as small spherical holes in the domain whose boundaries impose the activation onset. Under sufficient smoothness assumptions, an adjoint state can be defined through the shape derivative with respect to these boundaries, and therefore be exploited to optimize the EASs. This approach is efficient and can be applied to multiple pacing sites, although-in contrast to our approach-no topological changes are permitted by the formulation, impeding both a change in the number of EASs and movement from the interior to the boundary of the domain. Additionally, the coalescence of multiple sites needs special treatment in the shape derivative (not addressed in [22]). The viscous eikonal formulation, moreover, introduces a curvature-dependent conduction velocity, potentially strong at EASs, posing limitations on the radius of the spherical holes.
It is worth noting that all previous works, either based on the full [12] or a lumped [22,21] Purkinje network, require local measurements of activation times, e.g. endocardial maps, whereas GEASI can be applied to fit epicardial recordings and even the surface ECG. This aspect is relevant in view of noninvasive personalization of patient-specific models, as recently advocated [24]. Moreover, a major challenge in cardiac personalization addressed by GEASI is the estimation of the ground truth number of EASs. A possible solution is to consider a large number of EASs densely covering the earliest activation region, and successively removing sites according to some predefined rule [10]. For instance, an optimization procedure could determine the optimal activation onset of all sites, and then remove those with a very late onset. In a previous study, we optimized the initiation times along with the anisotropic conductivity tensors by manually deriving the Fast Iterative Method [9]. The large number of EASs results in a highly ill-posed inverse problem and consequently requires further regularization to successively remove initiation sites. While providing good results with respect to the measured activation times [20,10], this approach heavily relies on initial choice of EASs and the selected regularization strategy.

Notation
We denote the identity matrix by I and the space of symmetric and positive definite matrices in R d by Sym d×d . Throughout this work, B ζ (x) and B ζ (x) refer to the open and closed ball with radius ζ > 0 around x. We denote byX the interior of the set X.
[A] ij denotes the entry in the i-th row and j-th column of a matrix A ∈ R n1×n2 . Furthermore, A B if A − B is positive definite for A, B ∈ Sym d×d , and we set N := {1, . . . , N }. The set of unit vectors in R d is denoted by S d−1 . The Dirac measure of a set S is referred to as δ [S] . Further, we use the notation C 0 (X, Y ) for the space of continuous functions mapping from X to Y endowed with the norm · C 0 (X,Y ) , and we denote by C k (X, Y ) the associated space of k-times continuously differentiable functions equipped with the norm · C k (X,Y ) . We use the symbol C k,α (X, Y ) for the Hölder space with exponent α and norm · C k,α (X,Y ) . Finally, we denote by L p (X, Y ) and W m,p (X, Y ) the p-Lebesgue space and the Sobolev space of m-times weakly differentiable and p-integrable functions, and we set H m (X, Y ) = W m,2 (X, Y ).

Structure of the Work
In Section 2, we successively introduce the eikonal equation and the objective functional, which are the buildings blocks of GEASI. Based on the introduced algorithm, we present in Section 3 the topological gradient as well as the ECG fitting problem. Then, we elaborate in Section 4 on efficient discretization schemes and implementation detail of the proposed method. In Section 5, we consider the problem of estimating the initiation sites for different models-primarily in-silico experiments, but also one in-vivo experiment. Further aspects of future work are addressed in Section 6.

GEASI
Next, we introduce the GEASI method, which encompasses the following ingredients. In Section 2.1, we review the anisotropic eikonal equation and its associated Hamilton-Jacobi formulation. Subsequently, in Section 2.2 we analyze a general objective function involving the solution of the anisotropic eikonal equation from a functional-analytical perspective. Section 2.3 deals with the gradient computation of the distance function, which is later exploited in the aforementioned objective functional. Finally, all introduced concepts are combined in Section 2.4 to define GEASI.

Eikonal equation
We consider the computational domain Ω ⊂ R d for d ≥ 2, which in most cases represents the myocardium. Further, let E be the subset of N pairs for a priori given T min < T max and fixed N . Throughout this work, E is a set of EASs, where N is the number of EASs, x i and t i are the location and timing of the i-th site, respectively. Let φ E : Ω → R be the unique solution of the anisotropic eikonal equation with prescribed values on E, which is commonly referred to as the activation map. That is, φ E (x) is the first arrival time at x ∈ Ω of the propagating action potential.
where D ∈ C 1 (Ω, Sym d×d ) describes the anisotropic conduction. In the model, the anisotropy arises from the fiber alignment inside the heart [8]. Recall that Sym d×d is defined as the set of positive definite and symmetric d × d-matrices, which gives rise to the definition of the norm p D := √ Dp · p for p ∈ R d . Note that the assumptions already guarantee that for all x ∈ Ω and finite bounds 0 < λ * ≤ λ * < ∞. It is well known that the eikonal equation (1) admits a unique viscosity solution according to the theory of Hamilton-Jacobi equations [7]. The Lipschitz continuous solution of the eikonal equation φ E ∈ C 0,1 (Ω) is of the form where δ(x, y) denotes the geodesic distance given the length functional Thus, the induced Riemannian metric for two vectors We note that the infimum γ in (3) is actually attained, and by the geodesic equation we can even deduce γ ∈ C 0,1 ([0, 1], Ω) (see e.g. [25]). Indeed, in the definition (3), we first note that p D −1 (x) ≤ λ −1 * p 2 for all p ∈ R d and x ∈ Ω. Then, for any segment [x, y] fully contained in Ω we have that δ(x, y) ≤ When N > 1, all pairs (x i , t i ), (x j , t j ) ∈ E must satisfy the subsequent compatibility condition in order to ensure the existence of a solution. This fact is a direct consequence of (2), because non-compatible data can not exist w.r.t. the eikonal equation.
For the purpose of this work, this condition is not too restrictive, since we aim at identifying EASs rather than enforcing them. Interestingly, the condition is also physiologically sound: if a stimulus at some location x j is applied too late, e.g., right after the passage of an activation front originating from x i , it should not trigger another propagation, because the tissue is already depolarized. In fact, under such circumstances the activation time t j at x j would be larger than the travel time from x i , that is t i + δ(x i , x j ), clearly violating (6). The Hamilton-Jacobi formulation is essential for computing perturbations of E, which is conducted in the following subsection.

Objective functional
The overall objective of this work is the minimization of a given functional J : C 0,1 (Ω) → R depending on the activation map φ E with respect to E, i.e.
For instance, the objective could describe the minimization of a mismatch (in the least-squares sense) between the simulated activation and the activation detected from epicardial, as well as endocardial mapping (see Section 5). The objective functional can also involve the activation map implicitly: In Section 5.3, we utilize the mismatch between the recorded and simulated 12-lead surface ECG as a metric for optimization.
In what follows, we prove the existence of minimizers for (7) for varying N . To this end, we define for is a bounded function of its arguments.
Proof. Using (2), we immediately see that The Lipschitz continuity of δ as well as the compactness of Ω × U N imply the statement.
We note that Rademacher's theorem ensures the differentiability of Φ N almost everywhere. Non-differentiability with respect to x occurs for instance at x = x i , but also in the presence of front collisions. An immediate consequence of this lemma is the following Theorem 1 (Existence). If J is uniformly continuous, then the problem (7) admits at least one minimum.
Proof. The previous lemma and the uniform continuity of J imply the existence of at least one minimum. Proof. The first claim immediately follows from the definition of m N and set inclusion arguments. To prove the second claim, we assume that m N +1 = m N for some N and m N +2 < m N +1 . However, the choice x N +1 = x N +2 and t N +1 = t N +2 results in a contradiction. Corollary 1. If N is bounded from above by N max , then min N ≤Nmax m N has at least one minimum. Remark 1.
1. From a practical point of view, this corollary ensures that by adding new EASs, we either improve the objective function or we keep the same level of accuracy. This is also seen in the experiments in Section 5, where coalescence of two or more sites is observed if introducing too many EASs. (7) is in general not unique as it depends on the choice of J and on the order of the EASs. In principle, by permuting EASs we obtain the same value of the minimum. In particular, this symmetry induces a periodic partition of the set U N . Each partition is associated with a specific choice of the order of the EASs. From a numerical point of view, this may constitute a problem for methods based on random sampling. For deterministic steepest descent algorithms, the problem is mitigated by the fact that we rarely cross the boundary between two partitions, e.g., by swapping points, unless the two points coincide.

The minimum in
3. In general, we cannot take N unbounded with no further hypotheses on J . Suppose for instance that J is minimized by φ(x) = c for some constant c ∈ R. Then, inf N ∈N m N attains no minimum. Indeed, we cannot represent a constant function with (2) if E is only countable. However, we can approximate the constant with arbitrary precision with a sufficiently large number N of EASs.

Exponential Map
In what follows, we compute the Riemannian exponential map to derive an expression for the variation of the distance function. In particular, we discuss the relation of the derivatives of Φ N and the geodesic path. We briefly recall fundamental concepts in Riemannian geometry. Given x ∈ Ω and a tangent vector v ∈ V for a sufficiently small neighborhood V around the origin of the tangent space at x, the exponential map Exp x . In other words, the logarithmic map of y ∈ Ω identifies the tangent vectorγ(0) of a geodesic path γ emanating from x and ending at y. Proposition 2 (Variation of the distance function). Let x, y ∈ X , where X ⊂Ω is sufficiently small such that all points inside are connected by unique geodesics. Then the variation of δ(x, y) with respect to y with w = Log x (y) reads as Proof. Suppose that γ is a geodesic with respect to the Riemannian metric in (5) realizing the distance δ(y, x), i.e. γ(0) = y, γ(1) = x and Here, 0 < t 1 < · · · < t k < 1 are possible discontinuities of the geodesic curve and denote the one-sided derivatives from the left and the right, respectively. The derivative of γ with respect to the second argument is denoted by ∂ 2 γ. Since γ is assumed to be geodesic and smooth, the first two summands in (10) vanish.
In Proposition 2, we assumed uniqueness and smoothness of the geodesic curve, which is in general not ensured. In practice, the influence of geodesics violating these assumptions is negligible and thus in GEASI only consider (9) for all computations.
As before, let φ E be the solution of the eikonal equation with given This admits a natural definition of region of influences as Note that each point in the interior of R i is thus assigned to a single EAS through means of the geodesic distance. Furthermore, the derivatives of Φ N with respect to x i and t i at x ∈R i read as where we note that function is not differentiable on the boundary of the regions of interest. To compute the exponential map, we solve for each i = 1, . . . , N the following initial value problem for x ∈ Ω. The regularity and boundedness of D and φ {(xi,ti)} already imply the existence of solutions. Then, we define the piecewise geodesic path γ as is finite due to the assumptions regarding D for a small ζ > 0. The inclusion of the ζ-balls essentially circumvent problems related to the non-differentiability of γ i in the proximity of (x i , t i ). We note that by construction γ is a unit-speed geodesic for the length functional (4). In this case, we define the exponential map in the direction tγ(0) as Exp x (tγ(0)) = γ(t), and the logarithm Log x (γ(t)) = tγ(0) as its inverse. Note the the logarithm can efficiently be computed by tracking backward the geodesic from γ(t) to x. Remark 2. We note that points can belong to multiple regions of influence, at which the derivative of φ E might not be defined. However, due to the general functional-analytic setting the Lebesgue measure of these points is negligible.

GEASI Algorithm
In this section, we introduce the GEASI Algorithm to solve (7) using a gradientbased approach. Here, we restrict to the specific functional where Γ ⊂ Ω is a subdomain of Ω with a positive Lebesgue measure and φ ∈ L 2 (Γ, R) is a fixed square-integrable function. In numerical experiments, φ reflects the measurements on a known subdomain Γ, for which the quadratic mismatch on Γ between φ and φ with respect to the EASs is minimized. Examples of Γ include finite sets of points mimicking a contact recording map, the full endocardium/epicardium, or subregions of them. According to Sections 2.2 and 2.3, the optimization problem for N = 1 simply reads To employ a gradient-based approach, we see that following Proposition 2 the gradient of J with respect to x 1 simply reads as where γ x1→x (t) is the geodesic path from x 1 to x and r(x, is the residual. Optimizing multiple points simultaneously yields an average direction weighted by the residuals r on Γ. Figure 3 depicts how the velocity field shown in Figure 2 translates to a descent direction to optimize (15).  (15). Geodesics (white) originating from the single EAS x 1 to distinct points on Γ (left) and corresponding gradients computed with (9) (middle). The highlighted direction (red) coincides with the gradient in (16). Right: by iteratively applying a gradient-based scheme we determine the optimal (x 1 , t 1 ).
The extension to multiple EASs works similarly. A convenient formulation consists in splitting Γ into subdomains Γ i := R i ∩ Γ, each composed of those points activated by the EAS x i (note that the set of points belonging to multiple regions Γ i has Lebesgue measure 0). Then, the objective function reads as follows Clearly, the optimization procedure for a single EAS readily translates to the case of multiple sites. We found that rather than a simple gradient descent scheme, a Gauss-Newton optimization proved beneficial to reduce the overall number of required optimization iterations, resulting in the following update rule . (18) Here, } are the solutions of the previous iteration. To overcome local minima of the optimization problem (7) caused by non-unique solutions (see Remark 1) we additionally use an over-relaxation [27,28] with fixed β a = 1 √ 2 . The resulting Algorithm 1 iteratively linearizes and solves the problem using the computed gradient from δ to match a given measured activation. We remark that the gradient properly reflects infinitesimal changes of activation times on R i for each x i , but it is not capable of accurately capturing higher order effects like the change of R i . The experiments showed that rather than directly using E (k+1) from (18) as the new solution, it is beneficial to take a step-size β s < 1 and compute the convex combination of old and new solution according to this step size. For all experiments, we used β s = 1 2 . For further details of the numerical realization we refer the reader to Section 4.

Extensions of GEASI
GEASI is a versatile optimization algorithm, which can be extended in several aspects. In this section, we focus on two such possible extensions. First, the topological gradient estimation allows for an accurate estimation of the number of EASs. Second, we modify the original objective function of GEASI to fit a given ECG.

Variable number of EASs: Topological Gradient
So far, we assumed the number of EASs N to be fixed. Since the optimal number of EASs is in general unknown, we subsequently propose a method to approximate the optimal N . As a possible approach to estimate N (which is not conducted in this work) one could start with a large number of sites and successively remove distinct EASs that violate the constraint (6). However, this approach suffers from some major drawbacks: • several local minima can occur leading to a strong dependency on the initial guess, • enforcing (6) results in some numerical issues, e.g. dimension changes of the optimization problem and order of EAS removal.
In contrast, starting with a few (or even a single) EASs and subsequently introducing new EASs overcomes the above issues since according to Proposition 1 adding new sites does not increase the objective function. In what follows, we briefly recall the topological gradient, which is used to compute the infinitesimal expansion of splitting a single EAS. This expansion is exploited to estimate the decrease in objective function of adding a new site.
Consider the case of a single EAS, i.e. N = 1. The topological gradient is defined as the effect on the solution of the associated eikonal equation if splitting a single EAS x 1 into two new sites x 1 + εn and x 1 − εn in the direction of n ∈ S d−1 . We can directly infer from (2) that for E ε = {(x 1 + εn, t 1 ), (x 1 − εn, t 1 )}, where ε > 0 is sufficiently small. This topological operation divides the domain into two subdomains Ω − ε := {x ∈ Ω : We can now expand φ ε with respect to ε as follows where Φ 1 was defined in (8) and we used (9). In this case, we call the quantity the topological gradient.
A visual example of the topological gradient is provided in Figure 4. It is worth noting that the activation φ Eε continuously depends on the splitting distance ε. Hence, the topological operation of splitting an EAS does not introduce any discontinuities in the objective function. Moreover, we note that adding new optimal sites always decreases the objective functional unless ∇ x1 δ(x 1 , x)·n = 0. Therefore, we shall define a criterion for adding a split. The decrease in objective function of splitting a single site can be estimated as follows: Moving a single EAS in the direction n alters the activation times δ(x 1 , x) as shown. In contrast, splitting in the same direction n is similar to simultaneously moving a source point in both directions, and keeping only the shorter geodesic (right).
where r(x, Likewise, the effect of moving a source point in direction n is given by The ratio ν M,ε ν S,ε has proven to be a robust score for adding new sites, which is verified in the numerical experiments. In particular, if the ratio is below a certain threshold, then a new EAS is introduced.

Optimization using the ECG
The electrocardiogram (ECG) is the observed signature of the electric activity of the heart, which is measured at selected locations on the chest. Being routinely acquired and non-invasive, the ECG is the ideal candidate for inferring cardiac activation in a clinical framework. Here, we will introduce a method to reconstruct the EASs directly from ECG measurements. To this end, we exploit the methods presented in [5,6] to efficiently compute the ECG from activation maps of the eikonal equation. Finally, the quadratic mismatch of the computed and measured ECG is minimized, which yields optimal EASs. From a modeling perspective, we denote by Ω T ⊂ R d \ Ω the whole body domain excluding the heart cavity. The heart-torso interface Γ H := Ω T ∩ Ω is the boundary between the active myocardium and the rest of the body (for instance, endocardium plus epicardium), whereas Σ := ∂Ω T \ Γ H is the chest, on which the aforementioned electrical signal is recorded. We denote by T ⊂ R the considered time interval. In Figure 5, we outline how this setup would look in actual clinical measurements on the example of Lead I.
An equation for the torso potential can be derived from bidomain theory for the cardiac tissue and the balance of currents in the body (see e.g. [29]). Here, we consider the so-called pseudo-bidomain [30] or forward-bidomain model [31], in Figure 5: Exemplary setup for the lead-field of Lead I appearing in many ECG recordings. Ω T (torso) encapsulates heart domain Ω (heart domain), the left/right arm electrodes on Ω T are marked as x 1 and x 2 . The normal n is pointing from the surface of the heart domain Ω into the torso domain Ω T . The corresponding Z l is computed using (27), which is subsequently employed to obtain the ECG V l of Lead I.
which the parabolic and elliptic part of the bidomain equation are decoupled and can be solved sequentially. In this way, the transmembrane potential, denoted by V m (x, t, E), is not affected by the extracellular and torso potentials and can therefore be approximated independently. The resulting system of equations reads as follows: where the following quantities occur: u e (·, ·, E) : Ω×T → R is the extracellular potential in the heart, parametrized through the set of EASs E, The normal vector n at x ∈ Γ H points outwards, i.e. from the heart surface towards the torso, and is the outer normal vector for x ∈ Σ. The points x ± associated with x ∈ Γ H are obtained by taking the limit x ± ε = x ± εn for ε → 0.
The well-posedness of (23) follows from standard arguments for elliptic PDEs (see [32]). However, some care is required for the discontinuity across the hearttorso interface Γ H . Indeed, in order to render (25) meaningful, we need at least u T to be continuous on Σ for every t ∈ T. Let Ω = Ω ∪ Ω T ∪ Γ H be the domain modeling the whole torso (including the heart), and Following [29], we assume that Proposition 3. Under the above assumptions, the weak formulation of (23) given by is well-defined. In particular, there exists a unique solution up to an additive constant (the reference potential).
The proof directly follows from the Lax-Milgram theorem [32] by noting that V m (·, t, E) ∈ W 2,p (Ω) and u(·, t, E) ∈ W 1,p ( Ω) for p > d and all t ∈ T. Let x e ∈ Σ, e = 1, . . . , N E , be N E electrodes placed on the chest. A singlelead ECG recording is the potential difference between two such electrodes or, more generally, a zero-sum linear combination of the recordings. For instance, Einthoven's lead I is the potential difference between left and right arm electrodes. More generally, given electric potentials u T at the electrodes, the standard ECG is a vector-valued function V : T → R L given by where L is the number of leads and A is a L × N E real matrix defining the lead system, e.g. the 12-lead ECG. Since each row of A sums to zero, the matrix is not full-rank. For instance, the standard 12-lead ECG corresponds to the choice L = 12 (3 Einthoven leads, 3 augmented limb leads and 6 precordial leads) and N E = 9 (3 limb electrodes and 6 precordial electrodes), in which case A has rank 8. We remark that Morrey's inequality guarantees u(·, t, E) ∈ C 0 ( Ω), hence validating (25). Solving (23) is numerically costly for the standard 12-lead ECG, since we only evaluate u T at selected locations. Note that the system must be solved for every t ∈ T. Thus, we adopt the following integral representation of (25): where Z l : Ω → R are the lead fields (or Green's functions) satisfying the adjoint problem An informal derivation of (26) follows from the application of the second Green's identity to (23). As for (23), the solution is defined up to a constant. For a rigorous derivation accounting for the discontinuity in G, we refer the reader to [29, pp. 152 ff.]. A key observation is that the lead fields do not depend on t and E, making (26) particularly attractive for parameter estimation.
Next, we assign the transmembrane potential V m accordingly to a fixed waveform U : R → R shifted by the activation time φ E as follows We write the parametrized waveform as: which is visualized in Figure 6.  Table 1. The continuous formulation allows for an analytical derivation in (30).
Furthermore, the conduction velocity tensor D in the anisotropic eikonal equation in (1) is linked to the electric conductivity as follows: where β is the surface-to-volume ratio and α is a rescaling factor either experimentally estimated or obtained by solving the monodomain equation in a cable propagation setup [6]. Note that in all conducted experiments we assumed an equal anisotropy ratio G i = λG e , from which D = α 2 β λ 1+λ G i follows. The equal anisotropy ratio assumption simplifies the numerical experiments, but is not necessary for GEASI. All parameters adopted in this study are provided in Table 1.
We emphasize that φ E ∈ C 0,1 (Ω) only implies V m (·, t, E) ∈ W 1,p (Ω) and not V m (·, t, E) ∈ W 2,p (Ω) as required for t ∈ T and E ∈ U N . However, the aforementioned theory is still valid in this case with some major modifications that are beyond the scope of this work. Again, we refer to [29] and the references therein for further details.
In what follows, we intend to compute the sensitivities of the ECG with respect to the parameter set E ∈ U N . In the problem, only the activation map φ E appearing in the definition of V m depends on the parameters in E. Note that the chain rule straightforwardly implies The use of the aforementioned smooth waveform allows for a continuous analytical derivative ∂U ∂ξ . Details on the derivation of the term ∇ E φ E were already given in Section 2.3. Then, the derivative ∇ E V l is computed from (26) and reads as Finally, in this model the set of EASs E is computed from the measured ECG V l : I → R as follows: which is solved using the Gauss-Newton algorithm in a similar fashion to Algorithm 1. In particular, the update of the set E reads as follows with the modified objective functional The numerical integration in (32) is realized using the trapezoidal rule.
Remark 3. There are several numerical issues related to the optimization: 1. The waveform (28) is a rough approximation of a physiological action potential modelled the electrophysiology of a cell. The function U and the scaling parameter α may be simultaneously approximated from a generic ionic model by solving a 1-D propagation in a (possibly very long) cable with uniform coefficients. Alternatively, it is possible to show that (U, α) solves a nonlinear eigenvalue problem involving the ionic model [3]. (28) is actually not suitable to model the repolarization of the heart which is responsible for the T-wave. The reason is that the polarity of the T-wave, in general in accordance with the polarity of the QRS complex, can only arise from a heterogeneity in the action potential. Such heterogeneity might be introduced here, but it would be hard to reproduce the smoothing effect due to diffusion currents. Finally, the eikonal model is not suitable for the repolarization because, opposed to the depolarization phase, the repolarization front is of the same order of the size of the domain, impeding a proper perturbation analysis. In this work, the repolarization time is φ E (x) + APD, hence it satisfies the same equation as φ E , but with a shifted time.

Equation
3. Equation (30) requires higher order derivatives of V m and subsequently φ E . While we computed the derivative ∇ E V m as previously discussed, the computation of ∇ x V m is numerically achieved on the reference element.
4. It is important to mention that the gradient computation for the minimization of (31) is usually much more costly compared to optimizing the problem in the eikonal formulation from (15), since the size of Γ is much smaller compared to Ω. However, to compute ∇ xi,ti J we need the activation times and their derivatives in Ω, which necessitates the computation of the geodesics from each point of our domain to the EAS x i . The computational complexity is significantly larger than the complexity for (15). Further strategies to reduce additional computational costs are presented in Section 6.2.

Discretization
In this section, we elaborate on the discretization aspects for Algorithm 1, which encompasses the steps: over-relaxation of E (k+1) , solving the eikonal equation, computation of the geodesics and update of E (k+1) .

Solving the Eikonal Equation
The discrete function space for the eikonal equation is the space of volumetric Lagrange P 1 -finite elements defined on triangular (d = 2) and tetrahedral (d = 3) meshes discretizing Ω, respectively. Moreover, the discrete measurements in Γ are degrees of freedom (DOFs) of the mesh. Typically, finite element solvers require the initiation sites to coincide with DOFs of the mesh. However, since the original problem (7) expresses x i as a continuous quantity, we identify the DOFs of the actual element containing the activation site. Then, these DOFs are added to the Dirichlet boundary Γ D with fixed activation times given by t i + x − x i D −1 (xi) for x ∈ Γ D due to structural assumptions regarding the P 1 -finite element space. For the rare case of two or more initiation sites residing in the same element, we use the properties of (2) and (6) to compute the activation times.
In all subsequent computations, we employ the FIM [33,9] to solve the eikonal equation to account for the anisotropy.

Computation of Geodesics
In this work, we employ Heun's method (second order explicit Runge-Kutta scheme) to solve (12), which proved to be stable and efficient in numerical experiments. Due to the convergence of the ODE system to a stable node x i we terminate the iteration if the 2 -norm of two consecutive iterations is below 10 −10 . In practice, the ODE system is solved independently on each region of interest R i incorporating the whole set E. Since the gradient of the eikonal solution for the chosen discretization is a P 0 -finite element function (i.e. piecewise constant), we advocate a standard L 2 -projection onto the P 1 -finite element space [34]. Note that this projection can be realized by solving a linear system involving the mass matrix in P 1 . Since the boundary of Ω is in general curved, we project the geodesics back onto ∂Ω if they are outside of the domain after each update.
As remarked in Section 2.3, the gradient of the eikonal solution is discontinuous around each x i . To enforce regular gradients at each x i after the L 2projection of the previous eikonal solution φ E , we recompute the points with vanishing gradient by the subsequent variational problem with Tikhonov regularization for c ∈ d and balancing parameter λ > 0 as follows: Here, Ψ = (ψ 1 , . . . , ψ d+1 ) and {y j } d+1 j=1 are the collections of P 1 -basis functions and degrees of freedom associated with the element containing x i , respectively. Figure 7 depicts the effect of this regularization on the solution around x i .
The gradient computation in (9) is also sensitive to the choice of the step sizes of (12), which we choose as 5 · 10 −2 h with h being the average element size. As already described in Section 2.3, we compute the geodesic direction not directly at x i , but rather in a small ζ-neighborhood with ζ = 0.5h as advocated in (13). Numerically, the convergence of geodesics to this neighborhood is not ensured and non-converged geodesics (rarely occurring) do not affect the optimization.

Update of E (k+1)
Next, we optimize (18), where we have to ensure E ∈ U N . The constraint x i ∈ Ω is mesh-dependant and allows for no general analytical solution, potentially limiting the available optimization implementations. To overcome this hurdle, we use a proximal point algorithm enforcing E ∈ U N . The integration is realized using an exact simplex quadrature rule. In detail, we first compute the Moreau envelope of (18) with respect to the metric induced by . This particular choice of the metric [27] allows for an explicit solution to (34) given the projection onto U N . This method is usually referred to as Iterative Shrinkage and Thresholding (ISTA, [35]). In summary, the iteration step of the proximal point algorithm reads as Thus, the resulting optimal set is Note that the convexity of the projection depends on the convexity of the domain. In practice, hardly any cardiac mesh is convex, but nevertheless the proposed method generated reliable results for sufficiently small step sizes. The gradient direction of the Moreau-envelope is In this case, the unconstrained problem is solved using L-BFGS [36].

Numerical Results
Next, we present numerical results for various methods discussed above, where we focus on four setups to test GEASI on theoretical and cardiac problems: 1 The square domain presented in Figure 2 with a periodic conduction velocity field, where the measurement domain Γ coincides with the boundary of the surface. The initial initiation sites were chosen randomly for optimization w.r.t. activation times. For the optimization w.r.t. the ECG, we moved the target EASs further apart and used perturbations of the target EAS positions as initializations.
2 On a simplified 2D left ventricle (LV)-slice geometry with a transmural fiber rotation. Fiber and transverse intracellular conductivities were set to achieve a conduction velocity of 0.6 and 0.4 mm/ms, respectively. The measurement domain Γ is the outer ring of the domain, i.e. an epicardial slice. The initial initiation sites were chosen randomly for all related experiments.
3 On a clinically sampled, endocardial electrical mapping, recorded during intrinsic rhythm in a patient candidate to cardiac resynchronization therapy (CRT) and affected by a left bundle branch block. Data acquisition and the construction of the patient-specific anatomical model has been described in previous studies [6,10]. The measurements of activation were projected onto a patient-specific LV heart geometry. We mapped fiber orientations into the model using the approach described in [37]. Fiber, transverse and cross conduction velocities were set to 0.6, 0.4 and 0.2 mm/ms, respectively. Again, the initial initiation sites were chosen randomly, the target initiation sites are unknown.
In all experiments involving the computation of the ECG and lead fields, we assumed that intra-and extra-cellular conductivities are proportional to the tensor D, as explained in Section 3.2. Further numerical specifications of the aforementioned setups are listed in Table 2 and Figure 8.
In this work, we use a custom C++ implementation of the Fast Iterative Method [9] to solve the eikonal equation (1). Note however that the method is independent of the chosen eikonal solver and may benefit from higher order or smoother solutions of different solvers. A minimal working example for the method can be found on GitHub 1 , but is limited to the isotropic eikonal equation on structured grids using the Fast Marching Method [39,40]  1.5 · 10 4 10.7 · 8.9 · 9.5 1.9 2.5 4 1.08 · 10 5 10.3 · 8.1 · 12.6 0.66 3.5/18 (ECG) FEM. The ECG and geodesic computations, i.e. (26) and (12), and its Jacobian computation are calculated using the TensorFlow framework 2 , making use of available GPUs and enabling automatic derivation of V l with respect to V m . All computations were performed on a single desktop machine with an Intel Core i7-5820K CPU using 6 cores of each 3.30GHz, 32GB of working memory and a NVidia RTX 2080 GPU.

Activation Time Optimization
In Figure 9, we present the results of GEASI for the 2D experiments: In the first iterations of the square example, the EASs are moved to the center of the domain to promote a good overall fit during optimization. As the sites approach the center, fine details on the boundary can be fitted by minimizing the mismatch defining the optimal points. The idealized LV model additionally requires a non-convex projection since the fiber alignment favors movements on the endocardial wall. The optimization still works for this case, even though the problem in (34) becomes non-convex. Next, we concentrate on 3D experiments in Figures 10 and 11, for which we alter the number of initiation points for both models. Even though we can not ensure that the activation of the clinically acquired CRT patient can be described by the eikonal model with the simple rule-based fiber orientation, the results on the CRT masurements provide an overall low root-mean-square error (RMSE) between modelled and measured activation times. In the presence of a single EAS, the fit is (expectedly) sub-optimal since the activation requires a more complex activation pattern. With three or more EASs, we get a much better fit, evenly distributed throughout the ventricles, but additional initiation sites are moved from the septum to the LV. Note that there is no guarantee that the chosen rule-based fiber mapping can properly model the encountered ventricular activation. In such a case, the additionally employed initiation sites are able to compensate possible modelling inaccuracies. We can achieve even better results by sucessively increasing the number of initiation sites, but this only reveals the nature of the ill-posed problem: By increasing the complexity of our model, we can more closely approximate the presented activation map (cf. Proposition 1).
The trifascicular model has a higher resolution in comparison to the CRT   Figure 2 and exhibits no fiber orientation due to isotropy. Note that 3 was measured in-vivo and thus no ground truth is available.
model with an added fast conducting sub-endocardial layer, which is utilized in all longer geodesics from the measurements to the initiation sites. For this reason, it is important to properly project the geodesics in each iteration on the endocardium in a fast way. For further details of the actual implementation we refer the reader to Section 6.2. As a result, when using less EASs than in the ground truth we already achieve convincing numerical results, which is visualized in Figure 11. If we incorporate 6 initiation sites, we get a very good fit, even though one of the activation sites is deactivated before convergence due to (6). The three septal points are jointly modelled by two EASs accounting for the deactivated point. Adding points beyond the given ones did not yield any improvement as they are deactivated by other points during the optimization (not shown).

Topological Gradient
We also tried to estimate the correct number of EASs by using topological gradients (see Section 3.1), where we analyzed all 2D setups and the CRT patient. In this case, a splitting can only occur if the ratio ν M,ε ν S,ε is below 10 −1 (2D)/2.5·10 −1 (3D) and the maximum Euclidean distance of the position of two consecutive iterates among all EASs is smaller than 10 −2 h (with h being the average element  The minimizer of (21) is chosen by evaluating 360 (2D)/5625 (3D) directions, which are evenly distributed on the hypersphere. We additionally ensure that the splitting direction is feasible (i.e. it does not point outside the domain) by projecting the samples onto the mesh. To avoid two coalescing EASs inside one element after a split, the points are moved apart by 2h from the original site.
In Figure 12, we collected the results for all 2D experiments using this method and plot the ratio ν M,ε ν S,ε over the iterations. The first EAS is moved towards the center of the ground truth EASs, and subsequently several splits occur that closely match the ground truth sites. A similar behavior can be observed in the idealized LV model.
For the CRT patient in Figure 13, neither the ground truth EASs nor the fiber distribution and velocities in Ω are known. In total, the algorithm introduced 8 splits (i.e. N = 9), of which 4 are deactivated during optimization since they violated (6). Only those final EASs are shown in the right plot of Figure 13. Moreover, we can see that three main clusters are identified, where one initiation cluster is located at the upper part of the anterior septum. The optimization in this region is further complicated by the very thin wall of the 3D mesh, which likely causes the high number of splits. We highlight that constant (in time) split ratios are caused by temporarily deactivated EASs violating (6). To conclude, we get a tremendous fit with the presented measurement points despite the aforementioned model assumptions. Moreover, the topological gradient could be successfully applied to all 2D models leading to the correct estimate for the number of EASs and also matching the correct sites. The corresponding results in 3D provide a very low overall RMSE on the measurements.

ECG
In what follows, we present numerical results for the ECG optimization for both 2D experiments as well as the trifascicular model in a simplified fashion as a proof-of-concept. The ECG requires an additional full torso domain Ω T and the computation of the lead fields. For all experiments in this paper, we embedded all three in-silico experiments (i.e. 1 , 2 , 3 ) into a non-equilateral cube-torso without any additional organs and an overall torso conductivity of 0.2 mS mm −1 . The size of the cube-torso is proportional to the bounding-box of Ω. The computed lead fields are shown in Figure 14. In all cases, we generated a noiseless target ECG from the reference model setup with parameters and initiation sites already presented in Section 5. We optimize our model with random initialization with respect to this target ECG. Note that we do not focus on the generated ECGs' absolute potentials, since this heavily depends on the actual torso setup. Instead, we rather focus on the overall morphology of the ECGs.
To compute the lead fields in (27), our cube-torso is sampled using a structured regular grid of 100 d equidistant points, and the problem is solved with a finite difference scheme, which is sufficiently accurate since the lead field is  Figure 11: Results for the trifascicular experiment along with the RMSE (in ms) shown above the experiments for both Γ and Ω. The color-coded spheres indicate the observed activation times. The white and green circles represent the optimized and target EAS, respectively. The overall RMSE activation error is very low if using the correct number of initiation sites (N = 6), but we already obtain a good fit with fewer sites. evaluated far away from the singularity [41]. The lead fields are computed prior to the optimization since they remain constant. The ECG signals for the 3D models are mean-filtered with a small kernel of size ≈ 2 ms to improve accuracy. Optimization solely based on the ECG is frequently very challenging. However, with a proper initialization (x i , t i ), good fits for the ECGs can be computed. Figure 15 shows the optimization paths, as well as initial, target and optimized ECGs using the modified GEASI algorithm presented in Section 3.2 for the 2D examples, which are computed in approximately 2.5 hours each. The two potentials are a result of the two axis-aligned lead-fields (see Figure 14).
In the numerical experiments, it turned out that that the overall step size β s has to be chosen smaller compared to the activation timing problem. The morphology of the initial ECG and the optimized ECG differ by a large margin, making the fitting non-trivial. As a result, in both the square domain and the idealized LV experiment we are able to closely match the actual sites from which the target ECG was generated ( Figure 15, second row).
The trifascicular model in Figure 16 is computationally demanding since in each iteration step a computation of all geodesics is required, i.e. we need to solve ≈ 10 5 ODEs per iteration (for further details we refer to Section 6.2). As each initial EAS is randomly chosen, the initial ECG significantly differs from the target. Note that the 3D cube torso exhibits three axis-aligned leads. Since lead-X and lead-Z have the most prominent peaks, they have the largest effect on the resulting L 2 -error. After the optimization, these peaks were fitted by the algorithm by shifting most of the initiation sites to the LV and one to the anterior wall and septal region. The added difficulty with an activation featuring that many EASs is also apparent from the computed paths (white lines) which strongly vary during optimization. After termination, 4 of 6 sites are close to the ground truth sites defining the target ECG.

Discussion
In the previous section, we have experimentally demonstrated the broad applicability of the proposed GEASI method for a variety of problems. Despite

Splits
Split-Threshold Figure 13: Results for the topological gradient extension on the CRT experiment with a visualization analogous to Figure 12.
the convincing results there are still some issues related to our approach and alternative approaches, which will be addressed in future work.

Eikonal Equation
In this work, we rely on the anisotropic eikonal equation, but other versions thereof are also applicable. More specifically, several eikonal frameworks to model physical and medical processes have been proposed over the last three decades, which can be derived from either the monodomain or the bidomain equation using a perturbation argument [29]. The most common equation inferred from a first-order approximation of the monodomain equation is the anisotropic eikonal equation (1). The eikonal model originating from the bidomain model is slightly different and is based on a Finsler-type metric [42]. Second-order approximations lead to the curvature-eikonal, diffusion-eikonal and viscous-eikonal equations. In the curvature-eikonal model [3], the front velocity is corrected by the curvature of the front in the metric induced by the conductivity tensor. In contrast, in the diffusion-eikonal equation [4] a diffusion term is added to the right-hand side of (1). Finally, in the viscous-eikonal model [22], a squared eikonal equation is considered, which is corrected by a diffusion term.
Higher-order approximations have also been proposed, but are rarely used in practice [43]. The effect of higher-order terms is more pronounced in front collisions, at the boundary of the domain and in narrow channels, e.g., in scarred tissue. In practice, however, deviations from the classical eikonal model are minimal and is therefore widely accepted for personalization of cardiac models.
The distinction between these models is however important from the point of view of the EASs and is often dictated by the numerical method rather than the physiology. In the standard anisotropic eikonal model (1), EASs can be single points, whereas in the curvature-eikonal equation EASs are required to have a strictly positive Lebesgue measure. Note that the conduction velocity of a spherical front with small radius is significantly slower in the presence of higher-order correction terms.

Lead-Field X Lead-Field Y Domain
Square Torso -Lead-Fields Lead-Field X Lead-Field Y Domain Square Torso -Lead-Fields Lead-Field X Lead-Field Y Lead-Field Z Domain Figure 14: Setup for the ECG experiments showing the torso domain Ω T . The heart domain Ω is indicated by black lines for the 2D experiments and gray silhouettes for 3D. The streamlines visualize the lead fields. Note that the lead field for axis Z (green) is only present in the 3D experiments.

Runtime
The majority of the computational time is spent for solving the geodesics in (12), performed in parallel on the GPU. We highlight that the number of geodesics is proportional to the size of Γ in the original version (Algorithm 1) and proportional to Ω in the modified version (Section 3.2). The computation of all geodesics in both cases is performed in parallel on a GPU and therefore scales well with the mesh size. The bulk of computational time inside the ODE solver is spent on the projection of each ODE solution back onto the mesh and nearest neighbor computation. For the nearest neighbor computation, we implemented a custom KD-Tree implementation (publicly available on GitHub 3 ). For the projection operator, we extract the surface of the mesh, prior to the computation using the truncated signed distance function from VTK 4 . The K-nearest neighbor elements of the current positions of the geodesics are then queried to calculate the analytical projection onto all reference elements. The projection to all nearest neighbors is the minimum distance projection onto the mesh Ω.
Note that the time of each individual ODE solution in (12) depends on the length of the associated geodesic. Since adding new points can only shorten   Figure 15: Results of the 2D ECG optimization. Top row: temporal change of the positions of the EASs along with the ground truth. Bottom row: initial, final and target ECG for fitting. geodesic lengths (compare (3) and Section 3.1), more EASs will result in faster convergence and reduced computation time.
Solving the eikonal equation in (1) as well as the Gauss-Newton optimization in (15) only requires a minor portion of the computational time. As already mentioned, the activation time optimization is much faster compared to the remaining computations. In total, the experiments were finished within about 100 iterations only taking approximately 30 minutes and 90 minutes for 2D and 3D experiments, respectively. The experiments for the topological gradient behaves similarly regarding computational time. In contrast, the 3D optimization in the ECG problem requires approximately 12 hours.
To further decrease runtime, several approaches are possible: A custom GPU implementation to solve (12) along with the projection could significantly speed-up the optimization. Additionally, we often witnessed a collapse of many geodesic paths, especially in the trifascicular model, making subsequent computations redundant. An adaptive sampling from the measurement domain Γ   Figure 16: Results of the ECG optimization on the trifascicular model. Top row: optimized positions of EASs (white circles) along with temporal changes over the iterations (white lines). The green circles represent the target position from which the ground truth was generated. Bottom row: initial, final and target ECG. combined with a proper upsampling technique could increase performance at the cost of precision.
To improve performance for the 3D ECG optimization, we analyzed the convergence of the ODEs. Figure 17 shows a probability density function (PDF) of convergence of the geodesics γ over the number of required iterations using the trifascicular model with a single initiation site in the septum. Convergence in this case is defined as the first time two subsequent ODE iterations of (12) have a change of less than 10 −10 , i.e. γ(t k+1 ) − γ(t k ) < 10 −10 . We see that many of the computed geodesics converge very quickly, while points with a high geodesic distance need significantly more iterations before convergence. Our vectorized/parallel implementation to solve (12) exploits this fact to only include non-converged geodesics.

ECG
The ECG results demonstrated that GEASI can be used to fit a given ECG. However, one main problem is getting stuck in local minima. While the L 2 -error is relatively low in these minima, the morphology of the optimized and the target ECG differ a lot. One of the main reasons for this problem could result from the usage of the L 2 -error, which is not robust to transformations of the time series, such as time shifts. Better error measures for this type of optimization include dynamic time warping [44] and the Wasserstein distance [45]. Finally, different optimization algorithms could further help to overcome this issue.

Conclusion and Future Work
This paper introduced the novel GEASI method to find the optimal source points of an eikonal model with a special focus on electrophysiological examples. We showed that GEASI can model complex eikonal activations, either from the activation times directly, or by fitting a given ECG. For our model examples, we were able to identify most of the ground truth EASs along with times, and in the case of the topological gradient also the number. We were even able to model CRT measured data with only a few source points.
So far, we only assumed fixed conduction velocities and fiber distributions for all cases. In future studies we intend to estimate these parameters using the same procedure with only minor necessary modifications to (15). We note that past studies [20,18] have already shown that for the optimization of conductivities, additional regularization is crucial to decrease model complexity.
GEASI is inherently connected to the anisotropic eikonal equation through the Hamilton-Jacobi formalism. Thus, a possible extension of GEASI to reactiondiffusion models, possibly with non-local diffusion terms [46] and more complex boundary effects [30,47], is not trivial and probably requires a hybrid reactiondiffusion-eikonal approach [5].
All extensions of GEASI such as topological gradient and ECG optimization hold much promise for future applications in clinical real-world examples. GEASI faces several computational hurdles, many of which we already tackled in this study. We hope to further improve and expand GEASI-both methodologically and computationally-to enlarge the applicability to a wide-range of problems. Several pathological scenarios nicely fit in the GEASI framework, and the proposed method could potentially greatly improve the identification of the site of origin of premature ventricular contractions and monomorphic ventricular tachycardia [48]. GEASI could also be applied to improve planning of therapeutic interventions such as cardiac rhythm management with optimal placement and number of pacing leads [2]. Therefore, we believe that GEASI has the potential to significantly advance and improve personalized health care in the future.