Inverse localization of earliest cardiac activation sites from activation maps based on the viscous Eikonal equation

In this study we propose a novel method for identifying the locations of earliest activation in the human left ventricle from activation maps measured at the epicardial surface. Electrical activation is modeled based on the viscous Eikonal equation. The sites of earliest activation are identified by solving a minimization problem. Arbitrary initial locations are assumed, which are then modified based on a shape derivative based perturbation field until a minimal mismatch between the computed and the given activation maps on the epicardial surface is achieved. The proposed method is tested in two numerical benchmarks, a generic 2D unit-square benchmark, and an anatomically accurate MRI-derived 3D human left ventricle benchmark to demonstrate potential utility in a clinical context. For unperturbed input data, our localization method is able to accurately reconstruct the earliest activation sites in both benchmarks with deviations of only a fraction of the used spatial discretization size. Further, with the quality of the input data reduced by spatial undersampling and addition of noise, we demonstrate that an accurate identification of the sites of earliest activation is still feasible.


Introduction
Computational models of cardiac function are increasingly considered as a clinical research tool with the perspective of being used, ultimately, as a diagnostic Karl Kunisch was supported in part by the ERC advanced Grant 668998 (OCLOC) under the EU H2020 research program. Aurel Neic and Gernot Plank were supported in part by the Grant SFB MOBIS (FWF F3210-N18) and the BioTechMed Projekt "ILearnHeart".
B Philip Trautmann philip.trautmann@uni-graz.at Extended author information available on the last page of the article modality. Independently of which functional aspects are being considered, a key driving mechanism of cardiac electro-mechano-fluidic function is the sequence of electrical activation. Owing to its pivotal role, computer models intended for clinical applications must be parameterized in a patient-specific manner to approximate the electrical activation sequence in a given patient's heart. Anatomical (Demoulin and Kulbertus 1972;Ono et al. 2009) as well as early experimental mapping studies (Durrer et al. 1970), using ex vivo human hearts provided evidence that electrical activation in the left ventricle (LV), i.e. the main pumping chamber that drives blood into the circulatory system, is initiated by the His-Purkinje system (Haissaguerre et al. 2016) at several specific sites of earliest activation (root points) which are located at the endocardial (inner) surface of the LV. In a first approximation it can be assumed that the healthy human LV is activated at these root points by a tri-fascicular conduction system (Rosenbaum et al. 1969) consisting of three major fascicles referred to as anterior, septal and posterior fascicle. Owing to the fast conduction properties of the Purkinje network tissue patches surrounding root points are activated fast enough so that their activation can be considered instantaneous. Size and location of these patches as well as the corresponding instants of their activation are key determinants shaping the activation sequence of the LV. Since the His-Purkinje system is highly variable in humans, there is significant interest in inverse methods for identifying these sites, ideally non-invasively.
In general, non-invasive electrocardiographic imaging attempts to reconstruct the spatio-temporal behavior of the electrical sources of the heart from electrocardiograms recorded from the body surface by solving the inverse problem of electrocardiography (Gulrajani et al. 1989). Solving this inverse problem is complicated by the non-uniqueness of the relation between myocardial sources and their signature outside the heart, recorded in the form of extracellular electrograms. The vast body of research found in the literature can be broadly categorized based on the regularization techniques used to rule out solutions that are unlikely on physiological grounds (Tikhonov and Arsenin 1977) and the model used for representing the cardiac sources, with the predominant source models being transmembrane voltage-based (He et al. 2003;Wang et al. 2010), extracellular-potential based (Rudy and Burnes 1999;Bear et al. 2018), and activation/recovery-based (van Dam et al. 2009;Erem et al. 2014;Han et al. 2015;Janssen et al. 2018). These models have their pros and cons in terms of verifiability with experimental data, the domains in which sources can be reconstructed-on epicardial and endocardial surfaces or transmurally throughout the myocardial wall-and their accuracy in pathological scenarios such as the presence of infarcts (Wang et al. 2013) or more complex non-physiological activation patterns such as arrhythmias (Rudy 2013). For a comprehensive overview of these aspects of ECG imaging we refer to the recent review of Cluitmans et al. (2018).
In this study we propose a novel method for identifying these sites of earliest activation from activation maps measured at the epicardial (outer) surface of the heart. Such maps can be obtained non-invasively from body surface potential maps within clinical routine using inverse mapping systems such as CardioInsight (Ramanathan et al. 2004). Epicardial activation maps depend not only on location and timing of initial activation sites, but also on the orthotropic conduction velocities within the LV wall. Therefore, in patient-specific applications, conduction velocity tensors have to be identified using fast forward computational models (Zettinig et al. 2014;Marchesseau et al. 2013a, b), or biophysically detailed models (Potse et al. 2014). The propagation of electrical wavefronts in the LV is modeled based on the viscous Eikonal equation which is able to represent activation sequences and takes into account the dependency of conduction velocity on wavefront curvature. Identification of sites of earliest activation is achieved by solving a minimization problem. Initially geometries are chosen which represent the activation sites. Then they are relocated based on a perturbation field until a minimal mismatch between the computed and the given activation maps at the epicardial surface is achieved. The perturbation field is designed to reduce the functional subject to minimization during the relocation process. The proposed method is tested in two numerical benchmarks, a generic 2D unit-square benchmark serving the sole purpose of theoretical analysis, and an anatomically accurate MRIderived 3D human LV benchmark to demonstrate potential utility in a clinical context. For unperturbed input data, our localization method is able to accurately reconstruct earliest activation sites in both benchmarks with deviations of only a fraction of the used spatial discretization size. With the quality of the input data reduced by spatial undersampling and addition of noise, we demonstrate that an accurate identification is still feasible.
From a mathematical point of view the described problem can be interpreted as an inverse problem involving a non-linear elliptic PDE. On the activation sites ω i , i = 1, . . . , N an electrical depolarization wave is initiated which travels through the heart Ω = U \ ∪ N i=1 ω i . This is modelled by a nonlinear elliptic PDE, given by a viscous Eikonal equation, see Colli Franzone et al. (1990). The solution of the viscous Eikonal equation quantifies the arrival times of wave fronts at points in the heart Ω or on its surface Γ O . Since the wave is initiated on ∪ N i=1 ω i the arrival time is zero on ∂ ∪ N i=1 ω i and thus the viscous Eikonal equation has zero Dirichlet boundary conditions on ∂ ∪ N i=1 ω i and Neumann boundary conditions on the rest of the boundary of Ω. Given measurements of the arrival times on the surface of the heart Γ O the positions of the activation sites ω i are searched for. This inverse problem can be formulated as a shape optimization problem, see Delfour and Zolesio (2011) or Sokołowski and Zolésio (1992), in which the positions of ω i is optimized such that the misfit between the measured data and the solution of the viscous Eikonal equation on Γ O is minimal. We assume that the shape and number of activation sites is known and stays constant during the optimization. Thus only the locations of the activation sites are changed during the optimization. For the derivation of the shape derivative of the shape functional we use a technique which does not require the shape differentiability of the geometry-to-state mapping, see Ito et al. (2008) and Laurain and Sturm (2016). In order to apply this technique we first prove the wellposedness of the state equation. It is a nonlinear elliptic PDE which can be transformed to a linear one using the Hopf-Cole transformation, see Capuzzo Dolcetta (2003). The proof of the continuous dependence of the state on the data requires non-standard techniques. Furthermore we prove the wellposedness of the linearized and adjoint state equation using the weak maximum principle. In order to compute the shape derivative the averaged adjoint technique from Laurain and Sturm (2016) is used. In this manner we arrive at domain-based representation of the shape derivative, in contrast to the more common boundary-based representation, see Sokołowski and Zolésio (1992). This simplifies the numerical implementation of the shape derivative in a finite element environment, since only domain integrals need to be calculated. For the calculation of the perturbation field which is the basis for changing the geometry of the activation sites in an iterative gradient based algorithm a linear elasticity problem is solved in which the shape derivative enters as righthand side. To give a brief account of the contents of the paper, in Sect. 2 after the statement of the model on which our approach is based, we give its mathematical analysis, involving primal, tangent, and adjoint equations, and the shape derivative. The use of this information for numerical realization is described in Sect. 3. Finally Sect. 4 contains the two benchmark examples alluded to above.

Problem statement
Let U ⊂ R d , with d = 2 or d = 3, be a bounded domain with C 2,1 boundary, representing the cardiac domain. Within U we introduce N subdomains ω i with C 2,1 boundaries ∂ω i , which represent the volumes of the earliest activation sites, also denoted as activation sources. The union of ω i is denoted by ω = ∪ N i=1 ω i and its boundary by Γ = ∪ N i=1 ∂ω i . As such, Γ is the surface from which activation spreads into our computational cardiac domain Ω := U \ω. We have ∂Ω = Γ ∪ ∂U , and thus Ω is a bounded domain with C 2,1 boundary. In particular it is connected, but due to the holes it is not simply connected. Furthermore Γ ⊂ ∂Ω is closed. We set Γ N = ∂U , and further introduce the observatory boundary Γ O ⊆ Γ N , which in our application is given by epicardium of the heart. We consider the following minimization problem: subject to the viscous Eikonal equation in the form for some non-negative function g, and with The function T (x) represents the activation time, while the epicardial activation input data is denoted by z(x) which is assumed to be an element of L ∞ (Γ O ). The matrix M(x) models the squared cardiac conduction velocity (see Sect. 3.5). It is assumed to be symmetric and uniformly elliptic, i.e. there exists a α > 0 such that For the rest of this work we use the notation M ≥ α. The vector n denotes the outer unit normal vector on Γ N . The use of Eikonal equations is well-established to approximate the excitation process in the myocardium. We refer, for instance, to Colli Franzone et al. (1990) where a careful singular perturbation technique analysis with respect to the thickness of the myocardial wall and the time taken by the excitation wave front to cross the heart wall is carried out on the basis of the bidomain equations to arrive at various forms of Eikonal equations (Colli Franzone et al. 1990, Section 5).
Problem (1) falls in the class of inverse shape problems. For the numerical solution of (1) we require the shape derivative of J with respect to Γ in order to use it in a gradient decent method. As prerequisite we need to prove well-posedness of the state equation (2) which arises as PDE constraint in (1), and we analyze the tangent and adjoint equations.

Well-posedness of the viscous Eikonal equation
In this section, we discuss the well-posedness of the equation for some functions f , g specified later. Using the transformation T (x) = −ε log(w(x) + 1) this problem can be transformed into which is linear in the unknown w. Let us introduce the spaces

Moreover we set
For p > 1 let p its conjugate exponent. We introduce the positive and negative part of f defined by f + := max(0, f ) and f − := max(0, − f ) as well as the embedding constant c p > 0 of the embedding w L 2 p (Ω) ≤ c p w V . Next we require the following assumptions on the regularity of the data: Lemma 1 For every (M, f , g) satisfying (Ai), (Aii) and (Aiii), there exists a unique solution w ∈ V of (4). Moreover the solution satisfies w ∈ W 1, p denote the continuous trace operator onto Γ N . Using the embedding V → L 2d/(d−2) (Ω) and (Aii), it is easy to see that the integral Ω f wv dx is well defined for every v ∈ V . Due to the mentioned properties of the trace operator τ N and (Aiii) we can conclude that the boundary integral Γ N gwv ds is well defined. Thus we can formulate the weak form of (4) as To argue existence of a solution of (5) we use the Lax-Milgram theorem.
To prove the required coercivity in V we estimate for any w ∈ V using (Aii) and (Aiii) Thus we obtain coercivity and the existence of a unique solution w to (5). Moreover there exists a constant C > 0 depending on α and ε such that Next we argue additional regularity of w. For this purpose we consider the terms involving f w and gw as known inhomogeneities with w ∈ V . We show that the functionals F 1 (v) := Ω f wv dx and F 2 (v) := Γ N gwv ds are elements of (W 1, p (Ω)) * with p ∈ (1, 2) for d = 2 and p ∈ [6/5, 3/2) for d = 3. First we consider F 1 . We recall the embedding W 1, p (Ω) → Lq (Ω) withq = dp /(d − p ) = dp/(dp−d − p) andq = dp/(d + p). We prove that f w ∈ Lq (Ω). Using Hölder's inequality with r = (d + p)/d resp. r = (d + p)/ p we obtain Next we consider F 2 . We recall from Adams and Fournier (2003, Theorem 5.22 Here the restriction p ≤ 6 is necessary. Then assumption (Aiii) implies the assertion. Finally we get Moreover and a vector field f 2 ∈ L p (Ω, R d ), see Adams and Fournier (2003, Theorem 3.8). Thus the results from Troianiello (1987, Theorem 3.16) imply that w ∈ W 1, p 0 (Ω ∩ Γ N ) holds and the existence of a constant C depending on ρ M , ε and α such that These results are applicable since the Dirichlet part Γ of ∂Ω is closed.
In order to proof even higher regularity of w we use the following assumptions:

Lemma 2 Let Assumptions (Bi), (Bii) and (Biii) be satisfied. Then the solution of (4)
satisfies w ∈ C 2,δ (Ω) with 0 < δ < 1 given according to the data. Moreover there exists a constant C > 0 depending continuously on α, ε, ρ M , ρ f and ρ g such that Proof Theorem 3.28 (ii) and 3.29 (ii) from Troianiello (1987) can be applied, since (4) can be written as and Γ N is of class C 2,1 . This gives us the stated regularity and the corresponding a priori estimate. Next we define w + := max(0, w) This implies w ≤ 0 inΩ. Next we introduce the variableŵ = −(w + 1) which satisfies the equation If the solutionŵ were constant, it has to be equal to −1. However, in this case we haveŵ = 0 in Ω, which is a contradiction. We define O := max x∈Ωŵ ∈ [− 1, 0], see above. First we assume O = 0. Then Theorem 3.27 in Troianiello (1987) is applicable which states that such a maximum cannot be achieved on Ω ∪ Γ N . This is a contradiction. Thus O ∈ [− 1, 0) andŵ ∈ [− 1, 0). This implies the assertion.
For the rest of this work we fix a g ∈ C 1,δ (Γ N ) with g ≥ 0, 0 < δ < 1 and where the range of p is defined in Lemma 1. We define the set Proposition 1 There exists a constantc ∈ (0, 1) such that Proof We shall employ a compactness argument. For this purpose we argue that It is easy to see that M ≥ α holds. There exists a subsequence (M n k , f n k ) ∞ k=1 such that f n k converges for almost every x ∈ Ω to f . Thus f satisfies f ≥ β. On another subsequence of this subsequence there holds For the coordinate f this is a consequence of the estimates We remark that Lemma 1 is applicable for (M, f ) ∈ B and thus for elements of The solution w exists according to Lemma 1. The function δw satisfies the equation Then similar arguments as in the proof of Lemma 1 yield a constant C > 0 depending on ε, α and ρ M such that where p is specified in Lemma 1. The expressions involving w are estimated in terms of ρ g , ε, α and K . Thus there holds where h : R 2 → R is a continuous function with h(0, 0) = 0. Now, let (M,f ) be an arbitrary element in B Y . By Lemma 2 there exists a constant to conclude the desired result.
With the help of Lemma 2 we are able to define T = −ε log(w + 1) and calculate

Moreover we have on the boundary
We are now prepared to state the existence theorem for the state equation (3).
where C T only depends on B Y .
Proof Since existence of T was argued above only the estimate has to be proven. We know T = −ε log(w + 1), ∇T = − ε w+1 ∇w and Thus there holds wherec is the constant from Proposition 1 and K i only depends on B Y . This implies (10).

Well-posedness of the tangent and adjoint equations
Let Associated to the linearization of (2) we define the bilinear form Proof Let T be the solution of the state equation for M and f andT forM andf . Let w be the solution of (4) for M, f andw forM andf . Due to Taylor expansion of 1/x atw(x) + 1 the partial derivative of the difference δT := T −T satisfies the equation where δw := w −w and η(x) lies between w(x) + 1 andw(x) + 1. Due to Proposition 1 we have Now estimate (9) for δw in the proof of Proposition 1 with p = 6 and Lemma 2 imply the assertion.
Proposition 2 Let r ∈ (2, ∞) and F ∈ W 1,r (Ω) * . Then the linearized state equation Proof First we observe the following estimate Then we have Next we discuss the dependence of A −1 on M and T . First we remark that T depends on M and f . Thus we prove that the mapping (M, f ) → A is continuous from B Y endowed with the topology of C 0,δ (Ω, R d 2 ) × L 6 (Ω) to L(V , V * ). Let (M, f ) and (M,f ) be elements of B Y and A resp.Ã the corresponding operators. Then we estimate Finally we apply Troianiello (1987, Theorem 3.16, (iv)) which implies that v ∈ W 1,r 0 (Ω ∪ Γ N ) and whereĈ depends on ε, α, ρ M and C T .
Definition 2 For F ∈ V * we call ϕ ∈ V a solution of the adjoint state equation if it satisfies the equation A * ϕ = F or equivalently Theorem 2 Let r ∈ (2, ∞) and F ∈ W 1,r (Ω) * . Then Eq. (12) has a unique solution Proof From the proof of Proposition 2 it follows that A : V → V * is continuous and bijective. Thus A * : V → V * is also continuous and bijective. In particular we have (A * ) −1 = (A −1 ) * . So the equation A * ϕ = F has a unique solution ϕ ∈ V for every F ∈ V * and Then we apply Troianiello (1987, Theorem 3.16, (iv)) which implies that ϕ ∈ W 1,r 0 (Ω ∪ Γ N ) and whereĈ depends on ε, α, ρ M and C T .
Let us note that the strong form corresponding to (12) is formally given by

Shape derivative of J
We follow the notation and strategy in Ito et al. (2008) and Laurain and Sturm (2016).
for t ∈ [0, τ ]. We introduce and define the non-linear form e : After transformation to the reference domain Ω, the perturbed state equation can be cast as with the relation between T t and T t given by T t = T t • F t . Next we discuss the differentiability of A(t) and ξ(t). We shall use the notation where M k stands for the kth column of M.
Lemma 4 There holds where ξ (0) = div(h), and Proof Let x ∈Ω be arbitrary. The function ξ(t; x) has the form and ξ(t; x) = 1 + tr(Dh(x))t − (det(Dh 1 (x)) + det(Dh 2 (x)) + det(Dh 3 (x)))t 2 + det(Dh(x))t 3 , d = 3 where Dh i are the principal minors of Dh. Thus we have Thus the first assertion is proven. Let us turn to the differentiability of t → A(t). Since The derivative can be conveniently computed by its action on any Now let x ∈Ω be arbitrary and let t be so small such that t Dh * C(Ω,R d 2 ) < 1. Then there holds

Utilizing the product rule on A(t) = ξ(t)B * (t)M(t)B(t) leads us to (16).
The formulas for ξ and A also provide the following result.
The perturbed cost functional can be written as subject to e(t, T t , v) = 0 for all v ∈ V . Next we characterize the shape derivative at Ω in direction h. For this purpose we define the Lagrange functional

L(t, T t , p) = j(t, T t ) + e(t, T t , p)
for some p ∈ V and t ∈ [0, τ ]. We shall follow Laurain and Sturm (2016) to show that where T t solves (15) and ϕ t solves the averaged adjoint equation At first we characterize the right hand side of (19). First we observe that Since T t and T 0 appear linearly in (20), the averaged adjoint equation amounts to where [T t ] = 1/2(T t + T 0 ) ∈ C 2 (Ω).

Proposition 4 The averaged adjoint equation has a unique solution ϕ
In order to justify (19) we need the following technical lemma. (15) and of (20) for t ∈ (0, τ ]. Then we have

Lemma 6 Further let T t and ϕ t be the solutions of
Proof The first result follows from Lemmas 3 and 5. Let ϕ t be the solution of the averaged adjoint state equation (21) for A(t), [T t ] = 1/2(T t + T 0 ) and z. We define This follows from the fact that ϕ 0 ∈ W 1,r 0 (Ω ∪ Γ N ) → C(Ω) and ∇T 0 ∈ C 1 (Ω, R d ). Moreover the functional v → Ω A(t)∇δT ∇vϕ 0 dx is also a functional in V * , since δT ∈ W 1,6 0 (Ω ∪ Γ N ). According to the proof of Theorem 2 there holds with C > 0 independent of t. Moreover due to Theorem 2 there exists a constant c 1 > 0 depending only on B Y such that ϕ 0 This finishes the proof using Lemma 5.
We introduce the outer product v ⊗ w = vw * for v, w ∈ R d and the inner product G : N = trace(G N * ) for G, N ∈ R d×d . Now we have all necessary ingredients to prove the main result of this subsection. (18) satisfies

Theorem 3 The shape derivative d J (Ω, Γ ) of J defined in
Proof We apply Theorem 2.1 from Laurain and Sturm (2016). Thus we need to prove that lim t↓0 The functional J only depends on t through T t . Thus we have Thus we can estimate in the following way: Then Lemmas 6 and 4 imply the assertion. In order to calculate we recall Lemma 4 and in particular (16). We obtain Next we give a more usable formula for the shape derivative. For convenience we suppress the superscript for T 0 and ϕ 0 in the following. In particular we have

Practical implementation
In this section we describe the practical implementation of an algorithm utilizing the shape derivative D J for the reconstruction of the locations of the activation sites. We assume that these sites have the form ω i = B r i (x i ) with radii r i and midpoints x i , i = 1, . . . , N . For these activation sites we reconstruct the midpoints x i .

The state and adjoint state equations
Since the state equation is of nonlinear elliptic type which in practically relevant situations is posed on domains with challenging geometry, we propose to solve it using linear finite elements and a Newton method. For convenience we recall the state equation as In order to set up a Newton method we need to calculate the derivative of e, in particular we have The Newton equation is well posed, see Proposition 2. For a given solution T of the state equation, the adjoint state equation in the variable ϕ ∈ V has the form This is a linear elliptic equation of convection-diffusion type, which we again solve by linear finite elements.

Domain perturbation
While the overall source localization algorithm requires only a displacement of the current source locations, we still calculate a vector field for the perturbation over the whole domainΩ. This vector field h is chosen as the solution of the vector valued elliptic equation where S i , i = 0, 1 are defined in (24) resp. (23). We remark that h is defined on U and not only on Ω. The last equation is solved using linear finite elements. We also note that and thus h is a decent direction for J . Since we are only interested in the shift of the midpoints x i of the balls ω i , we average h over ω i , i = i, . . . , N , in order to get a shift of the midpoints.

Finite element solver implementation
The domain Ω is discretized using tetrahedral elements and linear Ansatz functions {ψ i }. As such, there are three linear systems to be solved at least once in each iteration of the source localization loop:

The domain perturbation equation
where S 0 and S 1 are defined according to respectively (24) and (23).
The linear systems are assembled and manipulated using the PETSc (Balay et al. 2017) framework. All three linear system are solved using the Boomer (Henson and Yang 2002). Algebraic Multigrid preconditioner in combination with the GMRES solver provided by PETSc. The linear solver in the Newton method is configured with a relative residual error tolerance of 10 −4 , while all other solvers use an absolute residual error tolerance of 10 −8 . The detailed solver settings are listed in the appendix.

Source localization
The goal of the source localization algorithm is to identify the midpoints x i , i = 1, . . . , N of the sources {ω i } that minimize our functional J . Our shape calculus based on shape derivatives does not allow for splitting or the creation of activation sites. For this purpose one has to resort to topological derivatives. We propose the approach depicted in Algorithm 1. Required inputs are some starting locations {x 0 i }, a user-specified, mesh dependent step-length (usually 1-3 mesh edgelengths), a step-length scaling parameter θ and a backtracking scale α. The symbol · denotes the Euclidean norm. The algorithm starts by initializing T 0 and J 0 . Then, while the tolerance condition on J k is not met, in each iteration of the while-loop it computes solutions to (27) and (28), updates the source midpoint positions and finally computes a new state solution to (25). If necessary backtracking is employed, and the next iteration begins.
For complex geometries, the step-length needs to be chosen small enough in order to prevent the sources from being moved out of Ω. Note, that only realizes an upper bound on λ i h k i , but this quantity is not bounded from below. Choosing θ > 1 improves convergence speed, as the λ i are scaled up to counteract the reduction of h k . In the case of overshooting, oscillations are reduced by backtracking.
According to the problem statement, the sources {ω i } are not part of the computational domain Ω. In each iteration k, all points ofΩ are moved based on the perturbation field h k , in particular the current source surface Γ k = ∪ N i=1 ω k i is moved. In practice it is easier to solve also the state and adjoint equations on U = Ω ∪ω with ω = ∪ N i=1 ω k i and apply the Dirichlet boundary values on wholeω. Then we only translate the logical representation of ω and thus the discretization of U is not perturbed. This prevents the need for re-meshing and implicitly enables the merging of any ω i without requiring special algorithmic treatment. Once λ i h k i is smaller than the average FE mesh edge-length, local refinement would become necessary. This however, is not within the scope of this work.

Algorithm 1 The source localization algorithm.
Choose initial mid-points x 0 Solve the adjoint state equation (27) for ϕ k . Solve (28) for h k .

Model parameters
The tensor parameter M contains the squared cardiac conduction velocity. In the depth of the human LV wall, conduction velocity is orthotropic due to numerous factors, with the most important ones being the geometry of myocytes and the non-uniform distribution of conduction-mediating proteins and sodium channels. The fastest propagation velocity v f is observed along the prevailing long axis orientation of myocytes, often referred to as "fiber orientation" f . Excitation spread within a sheet and along direction s, which is orthogonal to f , occurs at a lower conduction velocity v s , and even slower in a sheet normal direction n = f × s, at a velocity v n . Both orthotropic velocities as well as the principal axes { f , s, n} vary in space. In general, v f > v s > v n holds where the ratios are assumed as v f : v s : v n ≈ 4: 2: 1 based on experimental studies (Caldwell et al. 2009). As such, M is defined as The 2D benchmark in Sect. 4.2 will feature constant fiber-and sheet-directions f = (1, 0) * and s = (0, 1) * with varying (v f , v s ), while the 3D human LV benchmark in Sect. 4.3 will have constant velocities v f = 0.6 m/s, v s = 0.4 m/s, v n = 0.2 m/s and heterogeneous vectors { f , s, n}, computed by a rule-based method (Bayer et al. 2012). Further, in the human LV benchmark M(x) is an element-wise function. This makes the computation of S 0 impractical. While it would be possible to change the representation of M, this has not been pursued, since the terms involving S 0 have only a small impact on the shape derivative, see the comparisons in Sect. 4.2.
The parameter ε is calibrated by comparing the macroscopic velocity of propagating wavefronts generated by the viscous Eikonal model with physiological measurements such as the observed temporal delay between endocardial activation and epicardial breakthrough. Depending on a given trajectory relative to the used fiber field, macroscopic velocities fall into the range of local conduction velocities encoded in M, which themselves are based on experimental measurements (Caldwell et al. 2009).

Evaluation benchmarks
Two numerical benchmarks, a 2D wedge benchmark and a 3D LV benchmark, will be used to evaluate the proposed algorithm's ability to identify activation sources based on input boundary data.

Evaluation criteria
In both benchmarks we measure both the convergence of the current source locations {x k i } to the exact source locations {x i }, and the reduction of the functional J defined in (1). Thus the following evaluation criteria are used:
In this example we consider the noise free case. Thus the observed data z is generated by solving the state equation for T and restricting T to Γ N . In Fig. 1 we observe that the distances between the exact midpoints x i and x k i reach values below 10 −3 , more precisely d 1 = 1.7 · 10 −4 and d 2 = 2.6 · 10 −4 , after 100 iterations. These distances correspond approximately to the mesh size. On the right of this figure we can note that J k /J 0 attains a value of about 10 −7 . Figure 2 shows the trajectories of the points x k i as the iteration proceeds. We can see that the midpoints x k i do not move in straight lines. We expect that this is caused by interaction between the two activation sites, Moreover we see that the trajectories of the points x k i (Fig. 2) correlate to the main directions of the perturbation field h k .
In order to study the influence of the different parts of D J k we introduce the quantities: We clearly see in Fig. 2 that S 13 and S 14 are the dominating summands in |D J k h k |.
Thus it is justified to omit the terms S 01 and S 02 in the following benchmark. We also carried out tests with different choices for the conductivity tensor M and found nearly identical behavior provided that M is given by (29) with orthonormal vectors f and s. If the choice for M violates the orthonormality condition for f and s, then the numerical results may depend on the directions determined by the spatial relation between the exact activation sites and the initial guess, and the directions given by f and s.

3D LV benchmark
The 3D LV benchmark serves to gauge the potential of the proposed method in an envisioned clinical application which is geared towards localizing earliest activation sites from epicardial activation maps. In line with early experimental mapping studies (Durrer et al. 1970) on ex vivo human hearts we assume that there are three discrete sites of earliest activation located at the endocardial surface of the LV. In anatomical terms, these sites are located higher towards the base of the LV on the anterior paraseptal wall, a central area at the septal endocardium, and a posterior paraseptal area. Therefore the conduction system activating the LV is referred to as "trifascicular" with the three fascicles being referred to as anterior fascicle x af , posterior fascicle x pf and septal fascicle x sf . Each of these fascicles can be considered as a patch of tissue composed of a tight network of Purkinje fibers which are electrically coupled to the LV myocardium through so-called Purkinje-Ventricular junctions (PVJs). Owing to the fast conduction properties of the Purkinje fibers in these patches a large number of PVJs are located which activate one after the other with very short delays. Thus, these patches appear to activate simultaneously and are considered a fascicle and not a large set of individual PVJs. Further, due the short delays and the close spatial vicinity, it is still highly challenging today, even with invasive mapping devices recording signals with electrodes located in the immediate vicinity of PVJs, to identify individual PVJs. As such we do not expect that the identification of individual PVJs using data recorded at the epicardial surface is feasible.
While the presence of three fascicles is widely accepted and there general location is assumed to be known, the inter-individual variability and their exact location, size and relative timing is significant. Based on these considerations we assume the activation map (either measured or precomputed) on the epicardial surface as given input data for the localization of the three LV fascicles which we deem a plausible and sufficiently accurate general representation of the actual activation sources. Moreover, we simplify by assuming size and timing of individual fascicles as given and focus only on the identification of their location.
The discretized model of a human LV forming the computational domain Ω consists of 47,938 vertices and 245,611 tetrahedra, with an average discretization size of ≈ 1.5 mm. The observable surface Γ O is formed by the epicardial surface of the LV. We refer to Fig. 4. The source surface is Γ = ∪ 3 i=1 ∂ω i with ω i := B r =3mm (x i ), i = 1, 2, 3. Further, based on numerical tests with varying activation sequences we chose ε = 80 ms to obtain appropriate macroscopic conduction velocities which fall into the range of local conduction velocities encoded in M (see Sect. 3.5).
Motivated by real-world applications, we want z ∈ L 2 (Γ O ) to correspond to some error-prone data defined on a lower spatial resolution than the computational resolution of Γ O . To accommodate for this, the data z are generated as follows: is computed. 2. T r is sub-sampled at a set of 106 uniformly spaced sample points 3. A zero-average uniform noise is added: t i ← t i + ξ i 1 /|Ω| Ω T r (x)dx, i = 1, . . . , 106 with ξ i ∈ [− ξ /2, ξ /2]. 4. The data z is interpolated from the perturbed samples {t i } using distance-weighted interpolation.
We compare the following cases for different data selections: (RI) Using the reference data as input: z := T r | Γ O .
(II) Using only interpolated input: z is generated as described above with ξ = 0. (PI) Using perturbed and interpolated input: z is generated as described above with ξ = 0.3.
Starting at the initial locations the source localization Algorithm 1 is applied. In order to find the best achievable results, the algorithm is configured to only stop if J cannot be further reduced. All three cases executed in approximately 250 seconds on 10 cores of a workstation PC with two Intel Xeon E5645 (2.40 GHz) CPUs. Figure 5 shows the two evaluation criteria-the summed distances to reference location and the relative reduction of J -over the iteration count. For (RI), the algorithm terminates after 29 iterations with a relative error minimum of 3 · 10 −4 . The A B Fig. 6 a The trajectories traveled by x k i for the (RI) and (PI) cases. The (RI) trajectory is colored in green, while the one of the (PI) case is colored in red. The mesh vertices inside ball ω i , used for the reference solution T r , are displayed in red, while those inside the initial search ball ω 0 i are displayed in blue. b The adjoint solution ϕ k for the (RI) case, for the iterations k = 0, 10, 20, 28, respectively from left to right (color figure online) highest final distance to reference location is d 1 = 1 mm, which is well below the average FE edge-length of 1.5 mm. The discrete representation of the reconstructed activation sites closely match the desired reference sites.
In the (II) case, the algorithm stops after 29 iterations. The minimal relative error is 5.4 · 10 −3 . The highest final distance is d 1 = 4.3 mm, which is significantly larger than in the (RI) case. This indicates that the low-resolution sampling of T r | Γ O lowers the quality of our reconstruction. Also, the interpolation induces noise which impairs the reconstruction quality.
For (PI), the algorithm terminates after 30 iterations with relative error 1.6 · 10 −2 . The d i are similar to the (II) case, although slightly higher, with the highest final distance d 1 = 4.4 mm. This further hints that the low-resolution sampling has a much greater effect on the source locations than the error due to noise.
For all three cases, the final displacements x k i − x k−1 i are smaller than 0.2 mm, and therefore only a fraction of the mesh edge-length of 1.5 mm. As such, some mesh manipulation (e.g. mesh refinement, mesh deformation) would be necessary in order to apply the source displacement on the state and adjoint state problems. Since the mesh is not adjusted in the presented paper, this leads to a stagnation of the algorithm.  Figure 6 visualizes the source localization process by displaying the trajectories of x k i and the adjoint solution ϕ k during the source localization process. By comparing figure parts A and B, we observe that the motion induced by the field h is oriented from negative to positive regions of ϕ k , similar to the 2D benchmark in Sect. 4.2. Further, we see the diminishing absolute values of ϕ k over the iteration count. The final locations in Fig. 6a show, that even the worst localization (PI) still offers a good approximation of the general source location, well inside the uncertainty bounds of clinical parameters. Moreover, we carried out numerical tests with varying anisotropy ratios, see Fig. 7. In the LV benchmark, the convergence trajectory of one source varied significantly between the three choices of M. Numerical tests with significantly higher anisotropy ratios indicate, that a higher FE mesh resolution is required, particularly in the case of large displacements orthogonal to f .

Discussion
This study presented analysis and implementation of an algorithm for identifying sites of earliest activation in the LV from epicardial activation maps. The algorithm is posed as an optimization problem, where initial activation sites are chosen first to be then iteratively perturbed in order to minimize the mismatch between computed activation times and the activation maps given at the epicardial surface. We demonstrated well-posedness of all sub-problems, namely the viscous Eikonal equation, the tangent and adjoint equations and the perturbed state equation and characterized the shape derivative.
The theoretical results were verified by solving two benchmark problems, a 2D unit-square benchmark and a 3D human LV benchmark. For unperturbed input data, the localization method was able to accurately reconstruct the sites of initial activation. The largest deviations observed were 2.6 · 10 −4 and 1 mm, respectively, for the 2D and 3D benchmark. This was significantly smaller than the respective spatial discretization sizes of 4 · 10 −3 and 1.5 mm used in 2D and 3D benchmark, respectively. To probe the robustness of the method, the 3D benchmark was repeated using input data of reduced quality, that is, epicardial activation were spatially under-sampled and noise was added. These benchmark results showed, that the identification of earliest activation sites was still feasible, yielding a sufficiently accurate approximation of the general locations, comparable or better than the accuracy achieved with clinically used invasive endocardial mapping systems (Gepstein et al. 1997).
Several topics suggest themselves as possible extensions of the present work. The shape gradient is already set up to allow for a more realistic representation of the activations sites than those considered in the numerical realizations of these first benchmarks. Also, it can be of interest to incorporate different activation times by introducing inhomogeneous Dirichlet boundary conditions with unknown forcing terms. To allow for additional accuracy of the reconstruction of the evolution of the activation regions local grid refinement can be considered in future algorithmic efforts. Further it can be an interesting task to carry out the asymptotic analysis for ε → 0.

Limitations
While the benchmarks in this study demonstrate that the identification of sites of earliest endocardial activation from epicardial activation maps is, in principle, feasible with the proposed method, with regard to practical applications a number of restrictions apply. Out method makes various tacit assumptions which may not always hold in practice. Fiber arrangements are assumed to be known, following largely the patterns observed experimentally in the healthy LV (Streeter et al. 1969). With current technology fiber arrangements cannot be measured in vivo with sufficient spatial resolution,but suitable technologies under development (Scott et al. 2018) promise to lift this restriction in the future. Further, conduction velocities along the principal tensor axes were also assumed homogeneously throughout the LV, as velocities cannot be determined accurately in vivo, the chosen values were based on experimental observations (Caldwell et al. 2009). These values and their ratios may deviate from the experimentally estimation of v f : v s : v n = 3: 2: 1, and they may not be constant throughout the myocardium. Identifying the velocity tensor fields is therefore an additional complexity which is a related research topic (Marchesseau et al. 2013a) that has not been addressed in this study. A further limitation is the assumption that three sites of earliest endocardial activation exist. While this is physiologically motivated based on the notion that three main fascicles initiate activation in the healthy human LV endocardium (Durrer et al. 1970), this may not always be the case, particularly not under pathological conditions such as a left bundle branch block where the electrical activation of the LV may follow a markedly different pattern.
T , ϕ ms h mm M mm 2 /ms 2 ε ms