Digital Twinning of Cardiac Electrophysiology Models From the Surface ECG: A Geodesic Backpropagation Approach

The eikonal equation has become an indispensable tool for modeling cardiac electrical activation accurately and efficiently. In principle, by matching clinically recorded and eikonal-based electrocardiograms (ECGs), it is possible to build patient-specific models of cardiac electrophysiology in a purely non-invasive manner. Nonetheless, the fitting procedure remains a challenging task. The present study introduces a novel method, Geodesic-BP, to solve the inverse eikonal problem. Geodesic-BP is well-suited for GPU-accelerated machine learning frameworks, allowing us to optimize the parameters of the eikonal equation to reproduce a given ECG. We show that Geodesic-BP can reconstruct a simulated cardiac activation with high accuracy in a synthetic test case, even in the presence of modeling inaccuracies. Furthermore, we apply our algorithm to a publicly available dataset of a biventricular rabbit model, with promising results. Given the future shift towards personalized medicine, Geodesic-BP has the potential to help in future functionalizations of cardiac models meeting clinical time constraints while maintaining the physiological accuracy of state-of-the-art cardiac models.


Introduction
Cardiac digital twinning is a potential pillar of future's precision cardiology [1].A digital twin is a computational replica of the patient's heart, from its anatomy and structure up to its functional response.Commonly, it is obtained by fitting the parameters of a generic computational model to the clinical data of the patient.Here, we focus on cardiac electrophysiology, where we aim at personalizing a cardiac conduction model based on the eikonal equation, from non-invasive electrocardiographic (ECG) recordings [2,3].
Eikonal-based propagation models are an excellent basis for cardiac electrophysiology digital twins.They provide a good balance between computational efficiency and physiological accuracy for cardiac activation maps [4,5,6].When supplemented with the lead field approach, they can also provide accurate ECGs [7,3].Hence, the parameters of eikonal-based models could be fitted from non-invasive clinical data-surface potentials and cardiac imaging-in an efficient manner.This goal can be achieved using optimization methods based on stochastic sampling strategies [8,3], Bayesian optimization [9], gradient-descent [10], or deep learning approaches [11].
The numerical efficiency of the optimization procedure relies on two factors: the computational cost of the forward problem, and the number of samples required to achieve convergence within a given tolerance.Gradient-based methods can address the latter, if gradient computation is efficient and robust.For the former, many efficient eikonal solvers exist, but not all of them are suitable for automatic differentiation.A major criticality is that they internally solve a local optimization problem [12,13], one for each element of the grid and for possibly many iterations.This may represent a bottleneck in the computation of the gradient.A possible solution is to compute the gradient for the continuous problem, which can be shown is equivalent to the computation of geodesic paths [10,14], and then discretize.The optimize-thendiscretize strategy is however not ideal, because it may introduce numerical noise in the gradient and degrade the convergence rate.
In this paper, we introduce a novel inverse eikonal solver-named Geodesic-BP-and use it to efficiently solve the ECG-based digital twinning problem for ventricular models.To this end, we start by introducing a highlyparallelizable anisotropic eikonal solver for triangulations in arbitrary dimension.The key contribution is a novel approach for the solution of the local optimization problem and based on a linear convex non-smooth optimization method.Locally, the linear geodesic path within each simplicial element is approximated with a guaranteed bound, similar to many presented eikonal solvers [12].This solver can be implemented using current machine arXiv:2308.08410v2[math.OC] 5 Dec 2023

Figure 1:
The anisotropic eikonal equation is an efficient tool to compute the electrical activation of a heart (A), here shown on a public rabbit torso model.Many solvers utilize the relation of the eikonal equation to geodesic paths, which in the case of P1-meshes become piecewise-linear paths γ h between the initial conditions (xi, ti) and points on the surfaces x (B).By applying a non-linear transformation explained in Section 2.4, we can compute the measurable electrical activation on the torso surface and compute the mismatch against the measured potentials at sparse points (C).Implementing the approach in a differentiable manner using a backpropagation framework such as pytorch, we can efficiently calculate the gradient w.r.t. the initial conditions and use it to optimally choose (xi, ti) to match the measured potentials (D).learning (ML) libraries to efficiently calculate the solution of the eikonal equation.The forward solver is suitable for the computation of the gradient via the backpropagation algorithm, effectively enabling a discretize-then-optimize strategy (Geodesic-BP: Geodesic Back-Propagation).We will show that the backpropagation through such a solver is equivalent to piecewise-linear geodesic backtracking [15,16,10].Finally, we will demonstrate the potential of these methods for the problem of estimating initial eikonal conditions to match a given body surface potential map (BPSM) in a synthetic 2D study, and real 3D torso model.
In summary, the major contributions are as follows: 1.A novel, highly-parallelizable anisotropic eikonal solver, well-suited for GPGPU architectures and ML libraries.

A backpropagation-based gradient computation, in-
ternally related to the computation of geodesics.3.An efficient implementation of the method.4. We demonstrate how to use Geodesic-BP for the identification of cardiac activation sequence from the noninvasive ECG recordings. 5. We show that the inverse procedure is robust to model perturbations.
An outline of Geodesic-BP is provided in Figure 1.

Material and methods
In this section, the modeling assumptions for the electrical propagation inside the heart are introduced along with an efficient method for computing this activation in a parallelizable way.Then, we elaborate on the computation of the electrical activation given measurements on the torso leads, which ultimately leads to the formulation of the inverse ECG problem used throughout all experiments.

Cardiac activation
We model cardiac activation with the anisotropic eikonal equation (see [17,7,6]), which reads as follows where  [17].The onset of activation is determined by the set of initial conditions: , where x i ∈ Ω and t i are respectively the location and the onset timing of the i-th activation site.For convenience, we additionally define the norm in the metric D, i.e. ∥x∥ D = ⟨Dx, x⟩ for D ∈ S d .
Next, we discretize the domain Ω (assumed polytopal) by a simplicial triangulation.We denote by V = {v i } nv i=1 and E = {E j } ne j=1 respectively the set of vertices and elements of the triangulation.Each element E j ∈ E is a d-simplex (triangle in 2-D, tetrahedron in 3-D) with d + 1 vertices.We also define E j\i ⊂ ∂E j as the (d − 1)dimensional face of E j opposite to vertex v i .Note that E j\i is a (d − 1)-simplex (triangle in 3-D, segment in 2-D).For a general unstructured grid, we can approximate the activation time with a piecewise linear finite element function ϕ : Ω → R, see [18].(For the sake of simplicity, we use the same symbol as for the continuous case.)For the nodal values we use the shorthand notation ϕ i = ϕ(v i ).Following [12,19,20], we approximate D(x) = M −1 (x) as a piecewise-constant function and set D j = D(x)| Ej .
We numerically solve the anisotropic eikonal equation (1) by locally enforcing the following optimality condition for the nodal values: where ω i = {E ∈ E : v i is a vertex of E} is the set of elements containing v i and ϕ * j\i is the minimum of the following optimization problem: The optimality condition (2) is known as Hopf-Lax formula.It is possible to show that, under regularity assumptions on Ω, D, and the triangulation, the numerical solution of (2) correctly convergences to the unique viscosity solution of (1) when the triangulation is refined [21].As explained later, this principle is at the basis of Geodesic-BP.

Solution of the local problem
In what follows, we show how to efficiently and in a parallelizable way solve (3).We observe that eqs.( 2) and ( 3) consists of a series of constrained minimization problems, one for each element of the patch.Clearly, each local problem in ( 3) is convex, because the objective function is convex in x (note that ϕ is linear on E j\i ) and the constraint is a simplex, thus convex as well.Instead of relying on exact formulas, e.g., available in 2-D and 3-D from [20], we propose here a new dimension-independent algorithm for the solution of (3).With no loss of generality, we suppose that in (3) the vertices of E j are {v j1 , . . ., v j d+1 } and ordered such that the last vertex v j d+1 always corresponds to v i .Next, we consider the set of barycentric coordinates for the face E j\i : such that each point x in a face E j\i can be uniquely written as x = d n=1 α n v jn .Similarly, the P 1 function ϕ(x) on E j\i is given by the linear combination of the vertex values: ϕ(x) = d n=1 α n ϕ jn .Thus, Eq. ( 3) becomes where Φ j\i = (ϕ j1 , . . ., ϕ j d ) ⊤ is the vector of nodal values of ϕ| Ej , except for ϕ i , and A graphical representation of the local optimization problem is given in Figure 2. Eq. ( 5) now takes a common form of forward/backwardsplitting algorithms, where the function to minimize is a linear combination of two convex functions with one poten- Finding the minimum upwind origin x = n αnvj n can be cast into a non-smooth linear optimization problem, solvable using the FISTA algorithm (right).
tially being non-smooth.Such problems can be minimized using the FISTA algorithm [22], which alternatingly optimizes using a proximal step into the non-smooth function and a gradient step in the direction of the smooth function.The smooth part of ( 5) was thus chosen as h(α) = ∥Aα∥ 2 , and the non-smooth part as g(α) = ⟨α, Φ⟩+δ C d (α), where δ C is the indicator function on the set C such that The gradient of the smooth part h(α) is given by which is well-defined since h(α) > 0 for α ∈ C d and noncollapsed elements.The proximal step in the non-smooth part g(α) is then where L is a constant.The projection proj C d onto simplices is available in closed form [23, Figure 1].The ISTA algorithm consists of the following update step: An optimal choice of L in the algorithm is the Lipschitz constant of ∇h, which we can bound as where ∥•∥ in (8) refers to the spectral norm of the operator and h = 1 is a lower bound of h(x), for σ min (A) being the smallest singular value of A. Note that A is constant in each element, so the local Lipschitz constants can be efficiently pre-computed.
We summarize the solution of the local problem in Algorithm 1, where we also include the acceleration on Line 7 with β k = k−1 k+1 , as per the FISTA method, see [22].Note that Algorithm 1 works for arbitrary spatial dimension Algorithm 1: Solution of the local problem in Eq. ( 5) 1] Since the local problem is convex, it has a unique minimum α * .The rate of convergence is (see [22,Theorem 4.4]) Finally, we compute the global eikonal solution by iteratively applying Algorithm 1 to each element E j ∈ E and for all vertices v i ∈ E j , and then take minimum: see Algorithm 2.

Geodesic-BP
To utilize the computed eikonal solution in an inverse approach, we want to compute the gradients of the eikonal solution w.r.t.some quantity of interest, e.g., the initial conditions X 0 .More precisely, we are interested in the gradient ∇ (xi,ti) ϕ(x) and its efficient computation.This can be done analytically for (1), or algorithmically from Algorithm 2. The former (see e.g.[10]) relies on a optimizethen-discretize strategy.whereas the latter employs a discretize-then-optimize approach.Following the latter, we compute ∇ (xi,ti) ϕ(v) for all vertices v ∈ V and observe that ⊤ , (10) for x ∈ E j .Note that we neglect the derivative of α w.r.t.X 0 as a consequence of the first variation of the geodesic distance [24,Chapter 10].The nodal values ϕ jn in succession are then derived through (3), until we terminate at the simplex containing the initial condition x i ∈ E k .By defining x i in a continuous fashion (as was done in [10]), such that for all vertices v jn of the element E j containing x i it holds that the derivative readily follows: which can be effectively calculated using automatic differentiation of (11), or subsequent nodes (10) which are connected by a geodesic path to (x i , t i ).In particular, we use the back-propagation algorithm implemented in many ML libraries.More technical aspects are provided in Section 2.5.

ECG computation
We describe the torso with a domain Ω 0 ⊂ R d , whereas the heart is denoted by Ω, that is assumed to be well contained in Ω 0 .The full heart-torso domain is Ω T = Ω0 ∪ ΩT , and the torso surface is Σ = ∂Ω T .The electrodes are placed on the torso at fixed locations Xe = {x ei } N i=1 ⊂ Σ.Each lead V le (t) of the ECG is defined as the zero-sum of the electric potential at multiple locations.Thanks to the lead-field formula, we have the following representation: where and Z le : Ω T → R is the l e -th lead field, solution of the problem: The set X W contains the electrodes at the left/right arm and left foot, used to form the Wilson's Central Terminal (WCT).Eq. ( 13) can be solved once, since the problem is time-independent.The bulk conductivity of the torso, G i + G e , also depends on the extracellular conductivity tensor G e : Ω → S d .The derivation of this problem and exact formulation of (13) can be found in further detail in [5,25,10].
Next, we assume that the transmembrane potential is travelling wave with activation times dictated by ϕ(x): Since we are only interested in the depolarization sequence, we use the following template action potential, see [10,7]: where K 0 and K 1 are respectively the minimum and maximum transmembrane potential.Additionally, we precompute the discretized integral ECG operator (12) by means of a 2nd order simplicial Gaussian quadrature scheme.We used the library scikitfem [26] to assemble this operator for our P 1 mesh, automatically interpolating and integrating V m at the Gaussian quadrature points.This way, we can compute the leads as V le (t) = (BU (t − ϕ)) le where B : R nv×N is the mentioned discretized, precomputed lead field integration operator of (12), n v = |V | are the number of vertices and N is the number of leads.
For optimizing the inverse ECG problem, we want to fit V le of ( 12) against a measured Vle .We consider for that Algorithm 2: Geodesic-BP Compute all simplex upwind directions towards each vertex i inside element j 4 Compute for each element A j\i , L j\i #see Eqs. ( 5) and ( 8) With enabled gradient computations 5)

14
With enabled gradient computations where ϕ X0 is the eikonal solution for given initial conditions and N defines the number of measured/computed leads.Note that with B computed from (12), we can easily conclude for ξ = t − ϕ that ∇ (xi,ti) V le (t) = −B ⊤ ∂U (ξ) ∂ξ ∇ (xi,ti) ϕ, where ∇ (xi,ti) ϕ can be computed using Algorithm 2 as discussed in Section 2.3.Note that backpropagation through ∇ (xi,ti) V le will yield the equivalent result.The resulting gradients were then used for a gradient-based optimization utilizing ADAM [27].

Implementation aspects
There are a few practical remarks when implementing Algorithms 1 and 2 that are important to consider: When computing the local update in Algorithm 1, the number of iterations n f is the dominating factor for performance.Since the local problem is strictly convex with a precomputed Lipschitz constant L E , the exact error bounds (see (9)) towards the global minimum in (5) are known and could be used to compute n f on an element basis.Note that L j only depends on A j\i , which is strongly related to the mesh geometry and may increase as element quality decreases.
Furthermore, The upper bound presented in ( 8) is not tight and might lead to slow convergence of the local FISTA updates in Algorithm 1.In practice, we thus recomputed h by minimizing α * through Algorithm 1 with Φ = 0 and using the estimate h = h(α * )/1.1.
When considering Algorithm 2, note that the parallelization denoted in line 7 is more descriptive for CPU-parallel implementations.In ML frameworks, in contrast, the parallelization comes from the vectorization of single implemented operations, which might require re-ordering of some operations.Still, all operations inside the loop starting from line 7 can be parallelized, where only a special parallelization is required for Line 16 (scatter operation).This makes it well-suited for vectorized implementations on a GPU in frameworks such as PyTorch [28].
The optional performance optimization denoted in line 9 already significantly cuts down on the required memory for backpropagation, only computing updates when the solutions of simplex elements have changed.Checking for an ε change in activations before the update can also avoid unwanted gradient eliminations in cases where the old and new values of ϕ take the same value.
Note that the local problem in Algorithm 1 could also be solved using the exact formulation derived for lowdimensional cases (d ≤ 4) and fixed simplex dimensions as in previous works [12,20].In practice, however, the parallelization of this explicit solution is limited by the required condition to satisfy the constraint in (4).Such branching leads to unwanted performance losses in SIMD architectures when parallelizing over multiple elements.
Note that the number of activation sites is never modified by Algorithm 2. Some activation sites may however be overshadowed by the surrounding activation, thus effectively providing no contribution to the overall activation.When this happens, the gradient w.r.t. the overshadowed site is simply zero, thus the point is ignored by the algorithm.If an initial condition (x, t i ) is overshadowed as described at the final iteration of the optimization, we often refer to it as inactive.

Results
In this section, details of the dataset and results of applying our method for matching an ECG for both a synthetic 2-D and a real 3-D case are provided.

Simulation study
Setup To test and analyze our algorithm, we first start with a simple 2-D simulation study.We designed a simple setup consisting of an idealized torso with 8 electrodes distributed around the surface, with left and right arm electrodes composing Wilson's central terminal (WCT) (see Figure 3).Seven lead fields were then computed each w.r.t.WCT.

Figure 3:
The 2D simulation setup.We consider a simple symmetrical left ventricular slice Ω to be located in a heterogeneous torso ΩT , containing the blood pool ΩB and lungs ΩL, with different conductivities.8 electrodes xe i (blue circles) are circularly distributed around the torso, two of which are chosen as the WCT xW ∈ XW .
The LV, meshed at a resolution of 0.9 mm (≈ 2700 DoFs), contained 4 initial condition tuples (x i , t i ), located throughout half of the LV myocardium with timings between 2.5 ms and 10 ms.For computing the ground truth, we use another model with a 0.45 mm space resolution and perturbed conductivities in all sub-domains of the torso.As an initial guess, we consider 8 activation sites, all with 0 ms timing and evenly distributed at the transmural center.The optimization was performed using ADAM [27] for 400 epochs using a learning rate of 0.5, taking ≈ 5 min on a AMD Ryzen 7 3800X 8-core processor.

Results
The results of the optimization against the initial and ground-truth solution and parameters are shown in Figure 4. We can nearly perfectly fit the ECG and approximately recover the initial conditions.Overall, the total root mean squared error (RMSE) w.r.t. the eikonal solution ϕ is only 5.8 ms.We tested the algorithm also for a varying number of initial conditions and found that it is a hyperparameter of importance, but need not be chosen exactly: A cross-validation on the 2-D torso revealed that the worst match is achieved, by underestimating the number of initial conditions.In such a case, the achievable activations are not complex enough and the model is unable to faithfully replicate both the ECG and unobservable eikonal solution.Increasing the number of initial conditions beyond the ground truth's number is neither able to better match the ECG by a significant margin, but will also not overfit and significantly increase the error on the eikonal solution.The most common cause of high errors on the eikonal solution is related to getting trapped in local minima when optimizing for the ECG, circumventable by batched optimization (not utilized in any of the shown experiments).

Body Potential Surface Map Reconstruction
Setup The second anatomy is from a public rabbit ECGi model [29].The anesthetized rabbit was equipped with a 32-lead ECG vest during normal sinus rhythm and the signals were recorded with a frequency of 2048 Hz.The ECG was recorded for 10 seconds, filtered, and averaged over all beats to produce a single beat.From this beat, we extracted the QRS duration to match the biventricular activation.The biventricular mesh consisted of 508 975 DoFs and 2 522 501 tetrahedra, roughly equivalent to a mean edge length of 0.82 mm in a human heart (assuming a size factor of 3 between rabbit and human hearts).Note that we also tested the optimization with a lower resolution of 143 564 DoFs and 591 975 tetrahedra, roughly equivalent to a mean edge length of 1.25 mm.
We use the electrophysiological parameters from [30], including cardiac fiber orientation, lead fields Z l , and conduction velocity.In the conducted experiment, we chose K 0 = −85 mV, K 1 = 30 mV, τ 1 = 1 ms in (15).The extracted QRS sequence was chosen during the time interval t ∈ [100, 165] ms.We optimized for the initial conditions (x i , t i ) of the model, i.e.X 0 in Algorithm 2. For this purpose, we randomly distributed 100 points on the biventricular surface (x i ∈ ∂Ω) and with all timings initialized to the approximate Q onset t i = 20 ms.

Results
The experiment was run for 400 epochs using a learning rate of 0.5 on an Nvidia A100 GPU (40GB), taking approximately 1.8 h (or 0.7 h on the reduced resolution).The computed results can be found on Zenodo [31].
The results for optimizing (16) with 400 epochs of ADAM using Algorithm 2 can be seen in Figure 5.The intial and optimized activation are shown in the left upper and lower side of the image.is shown together with the comparison of all ECGs (right).The initially randomly sampled points evenly cover the domain Ω, to allow for a fast excitation of the whole domain, though this leads to overexpressed deflections in many leads initially.The optimization significantly changes the activation of the heart and is able to achieve a promising match of the ECG in many leads.The loss in (16) converges after ≈ 250 epochs, with only minor variations in the parameters thereafter.Of the initial 100 points, only 91 are still active at the end of the optimization, meaning that the corresponding t i is contained in the final solution and not earlier activated by closer initial conditions.

Discussion
We presented Geodesic-BP, an algorithm to efficiently solve the inverse ECG problem while being very efficient on a single GPU.We are able to fit the parameters of a patient-specific model only from the surface ECG and the segmented anatomy.The fitted model faithfully reproduced the recorded ECG in an animal experiment.In a synthetic case (Section 3.1), Geodesic-BP is able to recover the ground-truth solution with high accuracy, even in the presence of perturbation in the forward model and in the lead field.The resolution of the considered 3-D heart model consists of > 10 5 DoFs, which is sufficiently accurate for a high-fidelity eikonal solution in a human heart, as it would correspond to roughly 0.8 mm edge length.At full resolution, all computations only take 16 s.The considered cardiac parameter estimation problem is of high-relevance [32,33,1,34], yet only a few works are able to efficiently solve this problem using solely noninvasive ECG measurements.The considered heart model and ECG originates from a real-world rabbit rather than a human model since it is (to the knowledge of the authors) the only publicly available full torso and heart model, captured with this high-level detail and including preprocessed electrophysiological measurements.Such animal models share many properties of human cardiac electrophysiology [35].Furthermore, the algorithm and equations involved translate 1:1 to human cardiac electrophysiology.Likewise, the methods involved are not restricted to biventricular models, and can be extended to atrial, or whole heart models under sinus rhythm.We also foresee the future possibility of functionalizing publicly available human heart and torso models, currently missing parameterization [36].
The runtime in the full rabbit heart model of 1.8 h for the computed 400 iterations is still expensive but it reduces to 0.7 h for a coarser resolution (see Section 3.2).Two main advantages of our novel method are the high number of parameters that can be simultaneously optimized, al-lowing for many initial sites in contrast to other works [3,2], and having no spatial restrictions on the onset location [9].The first advantage is a consequence of utilizing backpropagation to compute the gradient w.r.t.our parameters, resulting in the runtime not being dependent on the number of parameters and the convergence being dictated by properties of the underlying function to optimize [37].Thus, compared to existing works based on derivative-free optimization with the same objective, our algorithm will scale much better in terms of spatially-varying parameters such as conduction velocities [38].The latter advantage of arbitrary onset locations follows from the use of the volumetric description of the eikonal equation and lead fields, from which we can compute parameters and the local activation times not only on the heart surface, but also transmurally.In conjunction, both of these advantages allow us to model potentially more complex activations than previously were able by other works.
From an implementation point of view, Algorithm 2 can be readily applied to optimizing the conduction velocity tensor D with little modification.Based on our previous research [32,16], however, optimizing conduction velocities either requires strong prior information in the form of regularization, or multiple activation sequences to create meaningful results.The question how well such information can be reconstructed from the BSPM, and the space of possible solutions w.r.t.initial conditions both remain topics for future research.
One major limitation of Geodesic-BP is that it cannot handle re-entry phenomena such as tachycardia and fibrillation, because it relies on the eikonal equation with focal boundary conditions.However, the eikonal model can be reformulated to allow re-entry, see [39] and it could be used in an optimization problem where the initial condition is to be identified from the ECG.Finally, we also foresee many practical improvements for reducing the computational time of the local solver as outlined Section 2.5.

Conclusions
This paper presents Geodesic-BP, an efficient method to compute piecewise linear geodesics, and demonstrated its performance for solving the inverse ECG problem in a limited parameter setting.By creating the link between backpropagation through a FEM-solver with geodesic backtracking, we were able to show that our previous FEMbased inverse eikonal models [15,16] are linked to our recent work [10] working with geodesic backtracking.Such inverse procedures are current research topics of highrelevance and efficient inverse eikonal models such as the present one may help in advancing patient-specific cardiac electrophysiology to practical real-world scenarios in the near future.

Figure 2 :
Figure 2: Conceptual visualization of a linear upwind discretization (5) on the example of a single tetrahedron (left).Finding the minimum upwind origin x = n αnvj n can be cast into a non-smooth linear optimization problem, solvable using the FISTA algorithm (right).

Figure 4 :FinalFigure 5 :
Figure 4: Result of the ECG optimization on the idealized torso experiment using 7 leads.The initial model (left) uses a simple simultaneous activation, evenly spread throughout the myocardium.The optimization reduces the number of initial conditions and is able to match both ϕ and V l closely (middle and right).The transparent red lines in the middle plot show the optimization path of initial conditions (xi, ti), which are inactive in the final solution.