A phase-field approach for the interface reconstruction in a nonlinear elliptic problem arising from cardiac electrophysiology

A phase-field approach for the interface reconstruction in a nonlinear elliptic problem arising from cardiac electrophysiology

Elena Beretta111Dipartimento di Matematica, Politecnico di Milano, P.za Leonardo da Vinci 32, 20133 Milano, Italy, {elena.beretta,luca.ratti,marco.verani}@polimi.it ,  Luca Ratti111Dipartimento di Matematica, Politecnico di Milano, P.za Leonardo da Vinci 32, 20133 Milano, Italy, {elena.beretta,luca.ratti,marco.verani}@polimi.it ,  Marco Verani111Dipartimento di Matematica, Politecnico di Milano, P.za Leonardo da Vinci 32, 20133 Milano, Italy, {elena.beretta,luca.ratti,marco.verani}@polimi.it
Abstract

In this work we tackle the reconstruction of discontinuous coefficients in a semilinear elliptic equation from the knowledge of the solution on the boundary of the domain, an inverse problem motivated by biological application in cardiac electrophysiology.

We formulate a constraint minimization problem involving a quadratic mismatch functional enhanced with a regularization term which penalizes the perimeter of the inclusion to be identified. We introduce a phase-field relaxation of the problem, replacing the perimeter term with a Ginzburg-Landau-type energy. We prove the -convergence of the relaxed functional to the original one (which implies the convergence of the minimizers), we compute the optimality conditions of the phase-field problem and define a reconstruction algorithm based on the use of the Frèchet derivative of the functional. After introducing a discrete version of the problem we implement an iterative algorithm and prove convergence properties. Several numerical results are reported, assessing the effectiveness and the robustness of the algorihtm in identifying arbitrarily-shaped inclusions.

Finally, we compare our approach to a shape derivative based technique, both from a theoretical point of view (computing the sharp interface limit of the optimality conditions) and from a numerical one.

\mdfdefinestyle

exampledefault rightline=false,bottomline=false, topline=false,leftline=false, backgroundcolor=white \mdfdefinestyleproblemdefault rightline=false,bottomline=false, topline=false,leftline=false, backgroundcolor=white \mdfdefinestylebadthingdefault rightline=false,bottomline=false, topline=false,leftline=false, backgroundcolor=white \renewbibmacroin: \DeclareFieldFormat[article]citetitle#1 \DeclareFieldFormat[article]title#1 \DeclareFieldFormat[article]pages#1 \DeclareFieldFormat[inproceedings]citetitle#1 \DeclareFieldFormat[inproceedings]title#1 \pdfstringdefDisableCommands \pdfstringdefDisableCommands

1 Introduction

We consider the following Neumann problem, defined over :

(1.1)

where is the indicator function of and

being and .

The boundary value problem (1.1) consists of a semilinear diffusion-reaction equation with discontinuous coefficients across the interface of an inclusion , in which the conducting properties are different from the background medium. Our goal is the determination of the inclusion from the knowledge of the value of on the boundary , i.e., given the measured data on the boundary , to find such that the corresponding solution of (1.1) satisfies

(1.2)

Since at the state of the art very few works tackle similar inverse problems in a nonlinear context, the reconstruction problem to which this work is devoted is particularly interesting from both an analytic and a numerical standpoint.

The direct problem can be related to a meaningful application arising in cardiac electrophysiology, up to several simplifications. In that context (see [book:sundes-lines], [book:pavarino]), the solution represents the electric transmembrane potential in the heart tissue, the coefficient is the tissue conductivity and the nonlinear reaction term encodes a ionic transmembrane current. An inclusion models the presence of an ischemia, which causes a substantial alteration in the conductivity properties of the tissue.

The objective of our work, in the long run, is the identification of ischemic regions through a set of measurements of the electric potential acquired on the surface of the myocardium. Indeed, a map of the potential on the boundary of internal heart cavities can be acquired by means of non-contact electrodes carried by a catheter inside a heart cavity; this is the procedure of the so-called intracardiac electrogram technique, which has become a possible (but invasive) inspection technique for patients showing symptoms of heart failure. We remark that our model is a simplified version of the more complex monodomain model (see e.g. [phd:tung], [book:sundes-lines]). The monodomain is a continuum model which describes the evolution of the transmembrane potential on the heart tissue according to the conservation law for currents and to a satisfying description of the ionic current, which entails the coupling with a system of ordinary differential equations for the concentration of chemical species. In this preliminary setting, we remove the coupling with the ionic model, adopt instead a phenomenological description of the ionic current, through the introduction of a cubic reaction term. Moreover, we consider the stationary case in presence of a source term which plays the role of the electrical stimulus.

Despite the simplifications, the problem we consider in this paper is a mathematical challenge itself. Indeed, here the difficulties include the nonlinearity of both the direct and the inverse problem, as well as the lack of measurements at disposal.

In fact, already the linear counterpart of the problem, obtained when the nonlinear reaction term is removed, is strictly related to the inverse conductivity problem, also called Calderón problem, which has been object of several studies in the last decades. Without additional hypotheses on the geometry of the inclusion, but only assuming a sufficient degree of regularity of the interface, uniqueness from knowledge of infinitely many measurements has been proved in [art:isakov_discontinuous] and logarithmic-type stability estimates have been derived in [art:ales]. Finitely many measurements are sufficient to determine uniquenely and in a stable (Lipschitz) way the inclusion introducing additional information either on the shape of the inclusion or on its size, e.g. when the inclusion belongs to a specific class of domains with prescribed shape, such as discs, polygons, spheres, cylinders, polyhedra (see [art:isa-po], [book:ammari-kang], [art:barcelo]) or when the volume of the inclusion is small compared to the volume of the domain (see [art:fried-vog], [art:cfmv]).

Several reconstruction algorithms has been developed for the solution of the inverse conductivity problem, and it is beyond the purposes of this introduction to provide an exhaustive overview on the topic. Under the assumption that the inclusion to be reconstructed is of small size, we mention the constant current projection algorithm in [art:amm-seo], the least-squares algorithm proposed in [art:cfmv], and the linear sampling method in [art:bhv] for similar problems. Although these algorithms have proved to be effective, they heavily rely on the linearity of the problem. On the contrary, it is possible to overcome the strict dependence on the linearity of the problems by aims of a variational approach, based on the constraint minimization of a quadratic misfit functional, as in [art:kv1987], [art:amstutz2005crack] and [art:ammari2012]. When dealing with the reconstruction of extended inclusions in the linear case, both direct and variational algorithms are available. Among the first ones, we mention [art:bhv] and [art:ikehata2000]; instead, from a variational standpoint, a shape-optimization approach to the minimization of the mismatch functional, with suitable regularization, is explored in [art:kalvk] [art:hett-run], [art:afraites] and [art:abfkl]. In [art:hintermuller2008] and [art:cmm], this approach is coupled with topology optimization; whereas the level set technique coupled with shape optimization technique has been applied in [art:santosa] and in [art:ito-kunish], [art:burger2003], [art:chan-tai], also including a Total Variation regularization of the functional. Recently, total-variation-based schemes have been employed to solve inverse problems: along this line we mention, among the others, the Levenberg-Marquardt and Landweber algorithms in [art:bb] and the augmented Lagrangian approach in [art:chen1999], [book:bartels, Chapter 10] and [art:bartels2012]. Finally, the phase field approach has been explored for the linear inverse conductivity problem e.g in [art:rondi2001] and recently in [art:deckelnick], but consists in a novelty for the non-linear problem considered in this paper.

Concerning the reconstruction algorithm for inverse problems dealing with non-linear PDEs, we recall some works related to sensitivity analysis for semilinear elliptic problems as [art:scheid], [art:amstutzNL], although in different contexts with respect to our application. We remark that the level-set method has been implemented for the reconstruction of extended inclusion in the nonlinear problem of cardiac electrophysiology (see [art:lyka] and [chavez2015]), by evaluating the sensitivity of the cost functional with respect to a selected set of parameters involved in the full discretization of the shape of the inclusion. In [art:BMR] the authors, taking advantage from the results obtained in [art:bcmp], proposed a reconstruction algorithm for the nonlinear problem (1.1) based on topological optimization, where a suitable quadratic functional is minimized to detect the position of small inclusions separated from the boundary. In [art:BCCMR] the results obtained in [art:BMR] and [art:bcmp] have been extended to the time-dependent monodomain equation under the same assumptions. Clearly, this type of assumptions on the unknown inclusions are quite restrictive particulary for the application we have in mind. In this paper we propose a reconstruction algorithm of conductivity inclusions of arbitrary shape and position by relying on the minimization of a suitable boundary misfit functional, enhanced with a perimeter penalization term, and, following the approach in [art:deckelnick], by introducing a relaxed functional obtained by using a suitable phase field approximation, where the discontinuity interface of the inclusion is replaced by a diffuse interface with small thickness expressed in terms of a positive relaxation parameter and the perimeter functional is replaced by the Ginzburg-Landau energy.

The outline of the paper is as follows: in Section 2 we introduce and motivate the Total Variation regularization for the optimization problem. In Section 3, after introducing the phase-field regularization of the problem and discussing its well-posedness, we show -convergence of the relaxed functional to the original one as the relaxation parameter approaches zero. We furthermore derive necessary optimality conditions associated to the relaxed problem, exploiting the Fréchet derivative of the functional. The computational approach proposed in Section 4 is based on a finite element approximation similarly to the one introduced in [art:deckelnick]. Despite the presence of the nonlinear term in the PDE it is possible to show that the the discretized solution converges to the solution of the phase field problem. We derive an iterative method which is shown to yield an energy decreasing sequence converging to a discrete critical point. The power of this approach is twofold: on one hand it allows to consider conductivity inclusions of arbitrary shape and position which is the case of interest for our application and on the other it leads to remarkable reconstructions as shown in the numerical experiments in Section 5. Finally, in Section 6 we compare our technique to the shape optimization approach: after showing the optimality conditions derived for the relaxed problem converge to the ones corresponding to the sharp interface one, we show numerical results obtained by applying both the algorithms on the same benchmark cases.

2 Minimization problem and its regularization

In this section, we give a rigorous formulation both of the direct and of the inverse problem in study. The analysis of the well-posedness of the direct problem is reported in details in the Appendix, and consists in an extension of the results previously obtained in [art:bcmp]. The well-posedness of the inverse problem is analysed in this section: in particular, we formulate an associated constraint minimization problem and investigate the stability of its solution under perturbation of the data, following an approach analogous to [book:engl, Chapter 10], but setting the entire analysis in a non-reflexive Banach space, which entails further complications. The strategy adopted to overcome the instability is the introduction of a Tikhonov regularization, and the properties of the regularized problem are reported and proved in details.

We formulate the problems (1.1) and (1.2) in terms of the indicator function of the inclusion, . We assume an a priori hypothesis on the inclusion, namely that it is a subset of of finite perimeter: hence, belongs to , i.e. the space of the functions for which the Total Variation is finite, being

endowed with the norm . In particular,

The weak formulation of the direct problem (1.1) in terms of reads: find in s.t., ,

(2.1)

being and . Define the solution map: for all , is the solution to problem (2.1) with indicator function ; the inverse problem consists in:

(2.2)

As it is proved in the Appendix, in Proposition 5, the solution map is well defined between the spaces and , thus for each there exists a unique solution .

We introduce the following constraint optimization problem:

(2.3)

It is well known that this problem is ill-posed, and in particular instability under the perturbation of the boundary data occurs.

A possible way to recover well-posedness for the minimization problem in is to introduce a Tikhonov regularization term in the functional to minimize, e.g. a penalization term for the perimeter of the inclusion. The regularized problem reads:

(2.4)

Then, we can prove the following properties:

  • for every there exists at least one solution to (2.4);

  • the solutions of (2.4) are stable w.r.t. perturbation of the data ;

  • if is a sequence of penalization parameters suitably chosen, then the sequence of the corresponding minimizers has a subsequence converging to a minimum-norm solution of (2.3).

Before proving the listed statements, it is necessary to formulate and prove a continuity result for the solution map with respect to the norm, which consists in an essential property for the following analysis and requires an accurate treatment due to the non-linearity of the direct problem.

Proposition 1.

Let satisfy the hypotheses in Proposition 6 or 7. If s.t. , then .

Proof.

Define ; then, subtracting the (2.1) evaluated in and the same one evaluated in , is the solution of:

(2.5)

where . Considering and taking advantage of the fact that and (by simple computation) , we can show that, via Cauchy-Schwarz inequality,

We remark that since . Moreover,

Thanks to Proposition 6, s.t. s.t.

where , and , which implies, thanks to the Poincarè inequality in Lemma 6.1,

Consider

since , then (up to a subsequence) pointwise almost everywhere, then also the integrand converges to . Moreover, , hence , and thanks to Lebesgue convergence theorem, we conclude that . Analogously, , and eventually , i.e. . Thanks to the trace inequality, we can assess that also . ∎

It is now possible to verify the expected properties of the regularized optimization problem.

Proposition 2.

For every there exists a solution of (2.4)

Proof.

Let be a minimizing sequence: then is bounded in and is bounded in (since is bounded and for all ). Thanks to the result of compactness for the space (see [book:ambrosiofuscopallara], Theorem 3.23), there exists a subsequence weakly converging to an element . Moreover, being weakly closed, . Since the weak convergence implies the convergence, thanks to Proposition 1 we can assess that in and in . Eventually, this proves that . Anologously, by semi-continuity of the total variation with respect to the weak convergence in BV, , and it is possible to conclude that

thus is a minimum of the functional. ∎

Even if the existence of the solution is ensured by the previous result, uniqueness cannot be guaranteed since the functional is neither linear nor convex (in general). We now investigate the stability of the minimizer of the regularized cost functional with respect to small perturbations of the boundary data. We point out that, due to the non-reflexivity of the Banach space , it is not possible to formulate a stability result with respect to the strong convergence; nevertheless, we can perform the analysis with respect to the intermediate convergence of functions. A sequence tends to in the sense of the intermediate convergence iff and .

Proposition 3.

Fix and consider a sequence such that in . Consider the sequence , where is a solution of (2.4) with datum . Then there exists a subsequence which converges to a minimizer of (2.4) with datum in the sense of the intermediate convergence.

Proof.

For every , we have that

Hence, and (and therefore ) are bounded, and there exists a subsequence such that both in and in . Thanks to the continuity of the map with respect to the convergence (in ) of and to the weak lower semi-continuity of the norm,

(2.6)

Hence, is a solution of problem (2.4). In order to prove that also , first consider that, according to (2.6),

hence

In addition, thanks to the continuity of , the first term in the sum admits a limit, i.e.:

which eventually implies that also . ∎

We finally state and prove the following result regarding asymptotic behaviour of the minimum of when .

Proposition 4.

Consider a sequence s.t. , and define the sequence of the solutions of (2.4) with the same datum but different weights . Suppose there exists (at least) one solution of the inverse problem (2.2). Then, admits a convergent subsequence with respect to the norm and the limit is a minimum-variation solution of the inverse problem, i.e. and s.t. .

Proof.

Let be a solution of the inverse problem. By definition of ,

Hence, is bounded, and since , is also bounded in and there exists a subsequence (still denoted as ) and s.t. . Moreover, , which implies that is a solution of the inverse problem (2.2), and

The lower semicontinuity of the norm with respect to the weak convergence, together with the continuity of the norm, implies that

for each solution of the inverse problem, which eventually implies that is a minimum-variation solution. ∎

Notice that, if the minimum-variation solution of problem (2.3) is unique, then the sequence converges to it.

The latter result can be improved by considering small perturbation of the data. By similar arguments as in proof of Proposition 4, one can prove the following

Proposition 5.

Let s.t. and let be such that and as . Suppose there exists at least one solution of the inverse problem (2.2). Then, every sequence , with , and solution of (2.4) corresponding to and , has a converging subsequence with respect to the norm. The limit of every convergent subsequence is a minimum-variation solution of the inverse problem.

Proof.

Consider a solution of the inverse problem. By definition of ,

(2.7)

In particular,

(2.8)

hence is bounded in and admits a subsequence (denoted by the same index ) such that : . Passing to the limit in (2.7) as ,

hence also

and by continuity of the solution map, we have that , which implies that is a solution of the inverse problem. By lower semi continuity of the BV norm (hence of the total variation) with respect the weak convergence and from inequality (2.8),

which allows to conclude that is also a minimum-variation solution of the inverse problem. ∎

Thanks to the results outlined in this section, it is possible to assess the stability of the regularized inverse problem:

In principle, one can obtain successive approximations of the solution of the inverse problem (2.2) by solving the minimization problem (2.4) with fixed . However, this approach would require to deal with several technical difficulties, namely the non-differentiability of the cost functional and the non-convexity of the space . This will be object of future research by adapting, e.g., the technique in [art:bartels2012] to the present context. However, in the sequel we follow a different strategy, namely introducing a phase-field relaxation of problem (2.4).

3 Relaxation

In this section, we formulate the phase-field relaxation of the optimization problem (2.4). Fixed a relaxation parameter , the Total Variation term in the expression of is replaced with a smooth approximation, known as Ginzburg-Landau energy or Modica-Mortola functional; moreover, the minimization is set in a space of more regular functions. We follow a similar strategy as in [art:deckelnick], with the additional difficulty of the non-linearity of the direct problem. In particular, we prove the main properties of the relaxed problem: existence of a solution is assessed in Proposition 1 and convergence to the (sharp) initial problem (2.4) as is proved in Proposition 3. Moreover, in Proposition 4, we describe the optimality conditions associated to the minimization problem, and compute the Frechét derivative of the relaxed cost functional, which is useful for the reconstruction purposes.

Consider and, for every , introduce the optimization problem:

(3.1)

The first theoretical result that is possible to prove guarantees the existence of a solution of the relaxed problem, employing classical techniques of Calculus of Variations.

Proposition 1.

For every fixed , the minimization problem (3.1) has a solution .

Proof.

Fix and consider a minimizing sequence (we omit the dependence of on ). By definition of the minimizing sequence, indepentently of , which directly implies that also is bounded. Moreover, being , a.e., thus and it is possible to conclude that is bounded in . Thanks to weak compactness of , there exist and a subsequence s.t. , hence . The strong convergence implies (up to a subsequence) pointwise convergence a.e., which allows to conclude (together with the dominated convergence theorem, since ) that

Moreover, by the lower semicontinuity of the norm with respect to the weak convergence, and by the compact embedding in ,

Moreover, using the continuity of the solution map with respect to the convergence, we can conclude that

Finally, by pointwise convergence, a.e., hence is a minimum of in . ∎

The asymptotic behaviour of the phase-field problem when is investigated in the next two propositions. First, we prove that the relaxed cost functional converges to in the sense of the -convergence, which will naturally entail the convergence of the corresponding minimizers. Before stating the next result, we have to introduce the space of the Lebesgue-measurable functions over endowed with the norm and consider the following extension of the cost funtionals: problem (2.4) is replaced by

(3.2)

whereas (3.1) is replaced by

(3.3)

It is now possible to formulate the convergence result, whose proof can be easily obtained by adapting the one of [art:deckelnick, Theorem 6.1].

Proposition 2.

Consider a sequence s.t. . Then, the functionals converge to in in the sense of the convergence.

Finally, from the compactness result in [baldo1990, Proposition 4.1] and applying the definition of -convergence, it is easy to prove the following convergence result for the solutions of (3.1).

Proposition 3.

Consider a sequence s.t. and let be the sequence of the respective minimizers of the functionals . Then, there exists a subsequence, still denoted as and a function such that in and is a solution of (2.4).

3.1 Optimality conditions

We can now provide an expression for the optimality condition associated with the minimization problem (3.1), which is formulated as a variational inequality involving the Fréchet derivative of .

Proposition 4.

Consider the solution map and let satisfy the hypotheses in Proposition 6 or 7: for every , the operators and are Fréchet-differentiable on and a minimizer of satisfies the variational inequality:

(3.4)

being

(3.5)

where and is the solution of the adjoint problem:

(3.6)
Proof.

First of all we need to prove that is Fréchet differentiable in : in particular, we claim that for it holds that , being the solution in of

(3.7)

namely, that

(3.8)

First we show that if , then . Indeed, the difference satisfies

(3.9)

with . Since and , and arguing as in the proof of Proposition 1, take in (3.9): we obtain

and again by Proposition 6

By (6.16) and by Sobolev inequality, eventually

hence .

Take now (3.9) and subtract (3.7). Define : it holds that

The second integral in the latter sum can be split as follows:

and in particular , where . Hence, chosen and exploiting again the Poincaré inequality in Lemma 6.1 and the Hölder inequality: