An Inverse Problem Involving a Viscous Eikonal Equation with Applications in Electrophysiology.

In this work we discuss the reconstruction of cardiac activation instants based on a viscous Eikonal equation from boundary observations. The problem is formulated as a least squares problem and solved by a projected version of the Levenberg-Marquardt method. Moreover, we analyze the well-posedness of the state equation and derive the gradient of the least squares functional with respect to the activation instants. In the numerical examples we also conduct an experiment in which the location of the activation sites and the activation instants are reconstructed jointly based on an adapted version of the shape gradient method from (J. Math. Biol. 79, 2033-2068, 2019). We are able to reconstruct the activation instants as well as the locations of the activations with high accuracy relative to the noise level.

Introduction

This work is concerned with an inverse problem in cardiac electrophysiology. In particular, the activation instants of the excitation wave in the myocardium are estimated from the arrival times of the wave at the epicardium. To briefly explain the problem we recall that the electro-physiologic activity of the heart is often modeled using the bidomain equations, whose numerical solution is very expensive. If one is only interested in the activation times T of the tissue, the bidomain model can be reduced to the simpler viscous Eikonal equation given, for instance, in the form Here the activation time T(x) is the time instance when the wave front passes through the point x in the domain Ω which describes the computational geometry of the heart. The epicardium of the heart is denoted by Γ and the boundaries of the activation regions (activation sites) by Γ. The matrix M describes the fiber orientation of the heart tissue and the values are the activation instants in the activation regions. On the basis of this model we formulate the inverse problem in the following form where z is the measured data on the epicardium. Problem (2) constitutes an inverse problem for the activation instants u. Similar inverse problems involving the Eikonal equation are considered in the context of seismic tomography, see for instance [23]. While in the analysis part we focus on reconstructing the activation instants from measurements of the activation time T on the surface of the computational domain Ω, in the numerical section we demonstrate that the activation instants and the location of the activation sites can be reconstructed simultaneously.

To briefly comment on the physiological background of this research, we point out that computational models of cardiac function are increasingly considered as a clinical research tool. For the understanding of the driving mechanism of cardiac electro-mechano-fluidic function, the sequence of electrical activations is of key importance. 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, which necessitates to solving inverse problems to identify patient specific parameters. Anatomical [5, 18] as well as early experimental mapping studies [6], 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 [11] 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 [21] consisting of three major fascicles referred to as anterior, septal and posterior fascicle. Size and location of these patches as well as the corresponding instants of their activation are key determinants shaping the activation sequence of the left ventricle. Since the His-Purkinje system is highly variable in humans, there is significant interest in inverse methods for identifying these sites and activation instants, ideally non-invasively.

As noted above it has become a standard procedure to rely for the mathematical description of the excitation process in the myocardium on various forms of Eikonal equations, see for instance [10, 17, 19]. These are reduced forms of the bidomain equations, a reaction diffusion system which describes the electrical activity of the heart. In [3, Section 5], see also [4], a singular perturbation technique 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 to arrive at various models for the Eikonal equation which differ by the nonlinear term. The two versions which are advocated in [3], and for which numerical comparisons are done there, are and . It is stated there that the model involving is better for wavefront propagation and collision. In the earlier work [15] we have used and solved the inverse shape problem of identifying the centers of the activation regions (spherical subdomains ω) from epicardial data z. An alternative approach for describing the wavefront propagation is based on the use of a function φ which at time t describes the level set of points for which φ(t, x) = 0 gives the position of the wavefront at that times, see [14]. We refer to [20] for an illuminating summary of wave propagation models involving Eikonal equations.

To briefly outline the paper, first we give a sufficient condition for the well-posedness of the elliptic PDE using the Schauder fixed point theorem and the maximum principle. The activation instants enter the state equation as constant Dirichlet boundary conditions on the surface of the activation regions. Then we calculate the gradient of the least squares cost functional with respect to these activation instants. It can be expressed in terms of the normal derivative of the solution to the adjoint state equation on the surface of activation sites. Therefore we also analyze the well-posedness of the adjoint and linearized state equations. Finally, we solve the least squares problem using the projected Levenberg–Marquardt method, see e.g. [9, 13] for a convergence analysis of this method.

In our first numerical experiments we consider only the reconstruction of the activation instants. In the second numerical example we perform the joint reconstruction of the activation sites and the activation instants. The activation sites are reconstructed by means of an adapted version of the shape gradient method introduced in [15] together with a projected gradient method for the reconstruction of the activation instants. The numerical examples illustrate the feasibility of the approach and are carried out on the 2D unit square with artificial data.

Problem Statement

Let , with d = 2 or d = 3 be a bounded domain and Γ = ∂U its boundary. In the physiological context it represents the cardiac domain. Within U we consider a family of open subdomains and we set Γ = ∂ω. These boundaries constitutes the surface from where the activation spreads. Then we define which is our mathematical and computational cardiac domain, with boundary . Note that Ω is connected but not simply connected. Let us choose a parameter ε > 0, and fix , which represents the epicardial input data.

With these specifications we consider the following problem: subject to the viscous Eikonal equation where β ∈ [0, 1], n is the unit normal on Γ, and . Further which is a bounded, closed, and convex set in . The function T stands for the activation time, and the matrix M models the cardiac conduction velocity.

Compared to the standard viscous Eikonal equation with the nonlinearity we introduce an additive parameter β > 0 in the model considered in (4). It is a mathematical tool which will allow us to guarantee that the control to state mapping u↦T is differentiable and that the adjoint equation is well-posed. For well-posedness of the optimal control problem in Section 3 the choice β = 0 is admissible. In Theorem 2 below we shall argue that the solution T of (4) as a function of β tends to the solution of (4) with β = 0.

Well Posedness of the Viscous Eikonal Equation and Existence of Optimal Controls

We assume that the boundaries of Ω are chosen such that the equation has a unique solution T ∈ H2(Ω) for any . Moreover we assume that and that M(x)v ⋅ v ≥ α|v|2 for a.e. x ∈Ω holds. Further, for any we assume the existence of g ∈ W2,6(Ω) with for i = 1, … , n, g vanishing in a neighbourhood of Γ, and , with c independent of u. For example can be chosen where the functions g are chosen as smooth bump functions which are equal to 1 on , vanish near Γ and have the property supp(g) ∩supp(g) = ∅ for all i, j = 1, … , n. Moreover for , we have for all i = 1, … , n and for all . In the subsequent developments (5) will be used with replaced by

Theorem 1

For ε > 0 sufficiently large (4) has a unique solution Moreover there exists a constant c, independent of , and β ∈ [0, 1] such that .

Proof

1. Existence: Let be fixed. Then we set with β ∈ [0, 1]. Since T ∈ V, g ∈ W2,6(Ω) and it follows that f(T) ∈ L2(Ω). Now let w ∈ H2(Ω) be the unique solution of with the estimate Thus we can define the operator , G: T↦w which satisfies the inequality with c(M, ε) independent of β ∈ [0, 1] and T. In the following we shall utilize Schaefer’s fixed point theorem in order to prove that G has a fixed point. At first we prove that is continuous and compact. Let {T} ⊂ V be a convergent sequence with limit T in V. We set w := G(T) and have according to (7). The compact embedding of H2(Ω) ∩ V in V implies the existence of a subsequence {w} and of a w ∈ V with in V. By taking the limit in the weak formulation of (6) we see that G(T) = w. Thus is continuous. A similar argument shows that is compact. In order to apply Schaefer fixed point theorem we have to further show that the set is bounded in V. Let T ∈ V be such that T = λG(T) for some 0 < λ ≤ 1. Then we have Multiplying this equation with T and integrating over Ω, we obtain by Young’s inequality and the fact that 0 < λ ≤ 1:

with c(M) independent of λ and ε. Thus if ε is sufficiently large, we have , for some constant independent of λ ∈ (0, 1] and β ∈ [0, 1]. Then Schaefer’s fixed point theorem can be applied to G and yields the existence of an element T ∈ V with G(T) = T which is a solution of (6). Setting T = T + g we have obtained a solution to (4), for which by (7) we have , with C(M, ε) independent of β ∈ [0, 1].

Moreover, since ∇T ∈ H1(Ω) and thus ∇T ∈ L6(Ω), and since also g ∈ W2,6(Ω) we have that and thus

2. Uniqueness: Let T ∈ W2,6(Ω), i = 1, 2 be two solutions of (4) and define δT = T1 − T2. Then δT satisfies the equation Let us define for the function It is easy to see, that holds. Indeed, in case the inequality is correct by the definition of B. Otherwise we have

Here we have used that defines a scalar product for the vectors . Alternatively we can note that B(x, v) is an element of the subdifferential of the convex function . Thus we have Consequently, where , since T2 is an element of W2,6(Ω). Then the maximum principle implies that δT ≤ 0 in Ω, see [25, Theorem 3.27]. Exchanging the roles of T1 and T2 in the above argument leads to δT ≥ 0 in Ω, and consequently to δT = 0, which implies the desired uniqueness. □

This proof is inspired from [8, Section 9.2, Theorem 5]. Henceforth it will be assumed that ε is large enough so that the solution to (4) according to Theorem 1 exists.

Theorem 2

We have where T denotes the solution to (4) as a function of β.

By Theorem 1 the family is bounded in H2(Ω) ∩ V and hence there exists a subsequence, denoted in the same manner, and such that in H2(Ω) and in V. Thus we can pass to the limit in to obtain that Moreover, by the trace theorem on Γ for i = 1, … , n. Now we set . By uniqueness, asserted in Theorem 1 we have , where T is the homogenous solution for β = 0 from Theorem 1, and thus the whole family converges to T0 in V. Moreover we have

for β → 0+. Since and defines an equivalent norm to the H2(Ω)-norm on H2(Ω) ∩ V, the claim follows. □

We close this section by asserting the existence of a solution to (3).

Theorem 3

There exists an optimal solution to problem (3).

Due to boundedness of U and boundedness from below of J, there exists a minimizing sequence {u} in U, satisfying Compactness and closedness of U imply the existence of a subsequence with a limit . Let denote the solutions to (4) with . By Theorem 1 the sequence is bounded in W2,6(Ω), and hence there exists a subsequence, denoted in the same manner, which converges weakly in W2,6(Ω) and, by Rellich’s compact embedding theorem, strongly in to some . We can now pass to the limit in (4) (with and ) to obtain that the pair satisfies (4). By (9) we have that . This concludes the proof. □

Well Posedness of the Linearized and Adjoint State Equation

For the practical realization of (3) the gradient of the cost is an essential tool for iterative solution schemes. For this purpose we investigate in this section the differentiability of the control to solution mapping, and we analyze the adjoint equation which will allow us to obtain a convenient expression for the gradient of the cost. These goals require us to assume β > 0, so that the nonlinear term in the state equation is smooth. An alternative procedure might involve proceeding as in [2, 16], at the expense of working with generalized differentiability concepts.

Throughout the rest of the theoretical part of this work T ∈ W2,6(Ω) with , and β ∈ (0, 1] are assumed. Further u, r and h are chosen arbitrarily in , L2(Ω) and , respectively.

We analyse the well-posedness of the following equations and

For this purpose we define the bilinear form by for any φ, v ∈ V. We recall the function g ∈ W2,6(Ω) defined in the previous section.

Definition 1

The function δT = v + g ∈ H1(Ω) is called a weak solution of (10) if v ∈ V solves the variational equation Analogously φ ∈ V is called a weak solution of (11) if it solves the variational equation

We introduce the operator and its adjoint by for all v, φ ∈ V.

Proposition 1

The operators and are isomorphisms. In particular there exists a constant C(M, T, ε) such that

The claims follow from a similar argumentation as in the proof of Proposition 2 in [15] using Garding’s inequality and the weak maximum principle. □

We introduce the space The space W is a closed subspace of H2(Ω), since the trace as well as the normal trace operator are continuous.

Proposition 2

Equation (10) has a unique weak solution which satisfies δT ∈ W and

First we define . We easily see that holds true. Thus Proposition 1 gives us the existence of v ∈ V satisfying (12) and we have the estimate

This implies that δT = v + g is the unique weak solution of (10). Moving the term L(T, v) to the right-hand side of (10), we conclude with standard elliptic regularity that δT ∈ W and that (13) holds. □

Proposition 3

Equation (11) has a unique weak solution which satisfies φ ∈ H2(Ω) ∩ V and

Proposition 1 implies the existence of a weak solution which satisfies the estimate Moving the div-term to the right-hand side of (11) and using which follows from the claim follows from standard elliptic regularity. □

Derivative of J

In this section we characterize the gradient of J using the linearized and adjoint state equations. It is an essential tool for every gradient based method, in particular for the numerical realization of problem (2), and will also be used for our numerical examples.

Lemma 1

There exists a constant C > 0 independent of β > 0 such that holds.

There holds Using the reverse triangle inequality for |⋅| we get □

Lemma 2

The function f : H2(Ω) → L2(Ω) defined by is Frechet differentiable with derivative for all β > 0.

By multiplication with the conjugate square root we get

utilizing Lemma 1 and . Then using the embedding H2(Ω)↪W1,4(Ω) and that we get □

Theorem 4

The operator , u↦T is Frechet differentiable and its derivative in direction is given by the solution δT ∈ W of (10) with u = δu for i = 1, … , N and r = 0.

We introduce the mapping defined by Using Lemma 2 it can be argued that E is Frechet differentiable. Moreover due to Proposition 2 the operator given by is an isomorphism. Let such that E(T0, u0) = 0. Then there exists a neighbourhood of T0 and of u0 and a Frechet differentiable implicit function S : U → V, u↦T with derivative given by δT = DE(T, U)− 1(0, δu). Since u0 is arbitrary, the result follows. □

Theorem 5

There holds where φ solves (11) for h = S(u) − z.

For each we have There holds where solves (10) with r = 0. □

A Projected Levenberg–Marquardt Method

We solve the inverse problem (3) based on a Levenberg–Marquardt strategy. Let be the orthogonal projection on U. In particular we iterate where 0 < λ ≤ 1 is the stepsize and d solves the problem The gradient of j is given by Thus we have to solve the equation Let be the matrix representation of .

Proposition 4

The matrix is positive definitive and there holds with w is the solution of (11) for h being the solution (10) for .

The formula follows from the exact same calculation as in the proof of Theorem 5, where we replace T − z by and φ by w. Moreover we have The corresponding equality implies on Γ. This fact, together with the unique continuation principle [1, 22] and uniqueness of solutions for the linearized state equation (10) imply that δu = 0. □

Numerical Example

In this section we present two numerical examples. In the first one we reconstruct the activation instants using the Levenberg–Marquardt method. As an alternative we could have used the iteratively regularized Gauß–Newton method. For a version with constraints see e.g. [12, 24]. In the second example we jointly reconstruct the positions of the activation regions and the activation instants using a combined shape gradient and projected gradient method.

Finding the Activation Instants

In this example, the computational domain U is given by the unit-square (0, 1) × (0, 1). We consider three activation sites ω = B0.1(x) whose midpoints are given by x1 = (0.5, 0.8)⊤, x2 = (0.2, 0.2)⊤ and x3 = (0.8, 0.4)⊤. Thus we have . The admissible set is given by . The observed data is given on the boundary Γ of U. While the domain is unrealistic from an applications point of view, the fact that the number of activation sites is known is physiologically reasonable. Indeed, early experimental work [21], using ex vivo human hearts provides evidence that electrical activation in the left ventricle is initiated at three discrete sites.

The domain U is discretized by 66049 vertices and 131072 triangles, which yields a discretization size of ≈ 4 ⋅ 10− 3. The state, linearized state and adjoint variable are approximated by P1 finite elements on the mentioned grid using the Fenics toolbox. The nonlinear system of equations for the state variable is solved using Newtons method. The linear equations for the linearized state and the adjoint state are solved using a direct method. Moreover we set ε = 0.1 and For the numerical computations we chose β = 0 and did not encounter any numerical instabilities due to term in the linearized state and adjoint state equations. In other examples instabilities might arise, which would force us to introduce β > 0. The exact activation instants are given by u† = (0, 0.1, 0.2)⊤. Then observed data z is generated by solving the state equation numerically for T for u†, restricting T to Γ and adding noise η. The used perturbance has the form where δ ≥ 1 and is a FEM-function with random coefficients on . The random coefficients are chosen from a standard normal distribution with zero mean and variance 1. The function is restricted to the boundary Γ and denoted by . Thus δ is the relative noise level. In this example we choose δ = 0.1 and δ = 10− 9.

After discretization of the state and adjoint state variable by linear finite elements, the projected Levenberg–Marquartdt Method can directly be applied. In every step of the iteration the matrix is calculated by solving the linearized state equation and the adjoint equation for all unit vectors resulting in a 3 × 3 matrix. Thus 6 linear PDEs must be solved numerically, see Proposition 4. Moreover the gradient of J has to be calculated by solving the nonlinear state equation and the adjoint state equation. So in total 8 PDEs have to be solved per iteration. The nonlinear state equation is solved by the Newton method. Since β = 0 the method can be interpreted as a semi smooth Newton method. The iteration is stopped by the discrepancy principle with τ > 1, see for instance [7]. In our experiments we choose τ = 1.1. The parameter α is set to 0.1.

In the case δ = 0.1 the discrepancy criterium is satisfied after 2 iterations with a final iterate , the state error and . In the case δ = 10− 9 the discrepancy criterium is satisfied after 5 iterations with a final iterate , the state error and . So we can observe that the activation instants are reconstructed very well relative to the noise level.

Finding the Activation Instants and Activation Regions

In this section we consider a similar scenario as before. But in addition to the activation instants we also reconstruct the position of the activation regions ω by determining the midpoints of ω. For this purpose we use the shape optimization approach introduced in [15] for the squared version of the Eikonal equation. Here we only modify the formulas developed in that work to fit our state equation. The shape derivative of J with respect to a smooth perturbation field h with compact support on is given by for any , where S, i = 0, 1 have the form

with the outer product v ⊗ w = vw⊤ for , the inner product G: N = trace(GN⊤) for , and where M stands for the k th column of M. Based on the shape derivative we calculate a perturbation field h by solving a linear elasticity equation of the form for γ > 0 and thus h is a decent direction for J. Since we are only interested in the shift of the midpoints x of ω, we average h over ω, i = i, … , N, in order to get a shift of the midpoints. The proposed method is of gradient type and thus we also update the activation instants based on the gradient calculated in Theorem 5. In particular we use a projected gradient method.

In the specific example we choose the exact activation sites as with , and . We denote . The exact activation instants are given by u† = (0, 0.1, 0.2).

We start the iteration at the initial points , and and initial times u0 = (0, 0, 0). Relative noise levels are chosen to be δ = 0.1, 0.01, 0.001 and the iteration is stopped using the discrepancy criterium.

In Table 1 and Fig. 1 we summarize our finding for the three noise levels δ. In particular we document the number of iterations K at which the discrepancy criterion is reached, the state error , the distance between reconstructed and exact positions denoted by and the reconstruction error . The reconstructed position of the three midpoints as well as the activation instants u are given for the respective noise levels by as well as

Table 1

Iterations and reconstruction errors for different δ

δ	Kδ	∥S(X_K_δ,u_K_δ)−z∥_L^2(Γ_N)	dK_δ1	dK_δ2	dK_δ3	|u_K_δ−u^†|
0.1	2	0.108	0.104	0.032	0.052	0.049
0.01	9	0.0104	0.021	0.006	0.013	0.013
0.001	52	0.001	0.003	0.001	0.008	0.005

Fig. 1

Top left: State error during the iteration; top right: Paths of the midpoints in U during the iteration; Bottom left: Distances between X and X† during the iteration; Bottom right: Error of the activation instants during the iteration

We conclude that the positions can be reconstructed with good quality relative to the noise level. Further tests in the noise free case showed that there is limit until which the state error can be reduced. This is caused by discretization effects.