On the inverse problem of detecting cardiac ischemias: theoretical analysis and numerical reconstruction

TL;DR

This paper develops theoretical analysis and numerical methods for detecting small inhomogeneities in myocardial tissue, aiding the identification of ischemic regions through boundary potential measurements.

Contribution

It extends asymptotic formulas for potential perturbations to 3D parabolic problems and introduces a topological gradient-based reconstruction algorithm.

Findings

Asymptotic formula for potential perturbations in 3D parabolic problems.

Numerical reconstruction algorithm demonstrates robustness.

Feasibility shown on idealized 3D heart geometry.

Abstract

In this paper we develop theoretical analysis and numerical reconstruction techniques for the solution of an inverse boundary value problem dealing with the nonlinear, time-dependent monodomain equation, which models the evolution of the electric potential in the myocardial tissue. The goal is the detection of a small inhomogeneity $\omega_\varepsilon$ (where the coefficients of the equation are altered) located inside a domain $\Omega$ starting from observations of the potential on the boundary $\partial \Omega$. Such a problem is related to the detection of myocardial ischemic regions, characterized by severely reduced blood perfusion and consequent lack of electric conductivity. In the first part of the paper we provide an asymptotic formula for electric potential perturbations caused by internal conductivity inhomogeneities of low volume fraction, extending the results published in [7] to the case of three-dimensional, parabolic problems. In the second part we implement a reconstruction procedure based on the topological gradient of a suitable cost functional. Numerical results obtained on an idealized three-dimensional left ventricle geometry for different measurement settings assess the feasibility and robustness of the algorithm.

Tables

Table 1. Table 1: Physical coefficients

$\nu$	$C_m$	$A^2$	$u_1$	$u_2$	$u_3$	$k_f^i$	$k_t^i$	$k_r^i$	$k_f^e$	$k_t^e$	$k_r^e$
$500\frac{m}{A}$	$0.1\frac{mAms}{c{m}^2}$	0.2	0	0.15	1	3	1	0.315	2	1.65	1.351

Table 2. (a) Results in presence of ischemia

Error	$N_p$	J	$mi{n}_{\Omega }G$
$0\%$	100	0.275	-0.5793
$1\%$	100	1.235	-0.589
$2\%$	100	4.128	-0.530
$5\%$	100	24.429	-0.589

Table 3. (a) Results in presence of ischemia

Error	$N_p$	J	$mi{n}_{\Omega }G$
$0\%$	100	0.275	-0.5793
$1\%$	100	1.235	-0.589
$2\%$	100	4.128	-0.530
$5\%$	100	24.429	-0.589

Table 4. (b) Results without ischemia

Error	$N_p$	J	$mi{n}_{\Omega }G$
$0\%$	100	0.000	0.000
$1\%$	100	0.964	-0.044
$2\%$	100	3.864	-0.105
$5\%$	100	24.148	-0.189

Equations

$$\{\begin{array}{ll}& \nu C_m{u}_t-\mathrm{div}(k_0\nabla u)+\nu f(u)=0,\mathrm{in}\Omega \times (0,T),\\ & \frac{\partial u}{\partial n}=0,\mathrm{on}\partial \Omega \times (0,T),\\ & u(0)=u_0,\mathrm{in}\Omega ,\end{array}$$

$$g(u,w)=-\beta u(u-\alpha )(u-1)-w,$$

$$\frac{\partial w}{\partial t}=\rho (u-\gamma w)\text{mbox}in(0,T),w(0)=w_0,$$

$$f(u)=A^2(u-u_1)(u-u_2)(u-u_3),u_i\in \mathbb{R},u_1<u_2<u_3,$$

$$\mathrm{dist}({\omega }_{\epsilon },\Omega \backslash K_0)\ge d_0>0.$$

$$|{\omega }_{\epsilon }|>0,\underset{\epsilon \to 0}{\lim }|{\omega }_{\epsilon }|=0.$$

$$\{\begin{array}{ll}& \nu C_m{u}_t^{\epsilon }-\mathrm{div}(k_{\epsilon }\nabla u^{\epsilon })+\nu {\chi }_{\Omega \backslash {\omega }_{\epsilon }}f(u^{\epsilon })=0,\mathrm{in}\Omega \times (0,T),\\ & \frac{\partial u^{\epsilon }}{\partial n}=0,\mathrm{on}\partial \Omega \times (0,T),\\ & u^{\epsilon }(0)=u_0,\mathrm{in}\Omega ,\end{array}$$

$$k_{\epsilon }=(k_0-k_1){\chi }_{\Omega \backslash {\omega }_{\epsilon }}+k_1=\{\begin{array}{ll}& k_0\mathrm{in}\Omega \backslash {\omega }_{\epsilon },\\ & k_1\mathrm{in}{\omega }_{\epsilon },\end{array}$$

$$\Omega \in C^{2+\alpha },\alpha \in (0,1),$$

$$u_0\in C^{2+\alpha }(\overline{\Omega }),u_1<u_0(x)<u_3\forall x\in \Omega ,\frac{\partial u_0(x)}{\partial n}=0\forall x\in \partial \Omega .$$

$$M_1:=\Vert f{\Vert }_{C([u_1,u_3])},M_2:=\Vert f^{\prime }{\Vert }_{C([u_1,u_3])}.$$

$$u_1\le u(x,t)\le u_3,(x,t)\in \overline{\Omega }\times [0,T],$$

$$\Vert u{\Vert }_{C^{2+\alpha ,1+\alpha /2}(\overline{\Omega }\times [0,T])}\le C,$$

$$u^{\epsilon }\in L^2(0,T;H^1(\Omega ))\cap C([0,T];L^2(\Omega )),u_t^{\epsilon }\in L^2(0,T;(H^1(\Omega ){)}^{\prime })+L^{4/3}(\Omega \times (0,T)).$$

$$\Vert u^{\epsilon }{\Vert }_{C^{\alpha ,\alpha /2}(\overline{\Omega }\times [0,T])}\le C,$$

$${\alpha }_1{u}^4-k\le f(u)u\le {\alpha }_2{u}^4+k,f^{\prime }(u)\ge -\lambda .$$

$${\int }_{\Omega }{u}_t^{\epsilon }v+{\int }_{\Omega }{k}_{\epsilon }\nabla u^{\epsilon }\cdot \nabla v+{\int }_{\Omega }{\chi }_{\Omega \backslash {\omega }_{\epsilon }}f(u^{\epsilon })v=0,\forall v\in H^1(\Omega ).$$

$${\int }_{\Omega }{u}_t^{\epsilon }v+{\int }_{\Omega }{k}_{\epsilon }\nabla u^{\epsilon }\cdot \nabla v+{\int }_{\Omega }{\chi }_{\Omega \backslash {\omega }_{\epsilon }}{u}^{\epsilon }v+{\int }_{\Omega }{\chi }_{\Omega \backslash {\omega }_{\epsilon }}\tilde{f}(u^{\epsilon })v=0,\forall v\in H^1(\Omega ).$$

$$\Vert z{\Vert }_{H^1(\Omega )}^2\le S(\Omega )\left(\Vert \nabla z{\Vert }_{L^2(\Omega )}^2+\Vert z{\Vert }_{L^2(\Omega \backslash {\omega }_{\epsilon })}^2\right),\forall z\in H^1(\Omega ),$$

$$a_{\epsilon }(u^{\epsilon },u^{\epsilon })={\int }_{\Omega }{k}_{\epsilon }|\nabla u^{\epsilon }{|}^2+{\int }_{\Omega \backslash {\omega }_{\epsilon }}(u^{\epsilon }{)}^2\ge S\Vert u^{\epsilon }{\Vert }_{H^1(\Omega )}^2,$$

$${\varphi }_n\in C^1(\overline{\Omega }),0\le {\varphi }_n(x)\le 1,\forall x\in \overline{\Omega },{\varphi }_n(x)=1,\forall x\in \overline{\Omega }\backslash {\omega }_{\epsilon },\mathrm{and}{\varphi }_n\to {\chi }_{\Omega \setminus {\omega }_{\epsilon }}\mathrm{in}{L}^{\infty }(\Omega ),$$

$$\{\begin{array}{ll}& u_t^n-\mathrm{div}(((k_0-k_1){\varphi }_n+k_1)\nabla u^n)+{\varphi }_{n}f(u^n)=0,\mathrm{in}\Omega \times (0,T),\\ & \frac{\partial u^n}{\partial n}=0,\mathrm{on}\partial \Omega \times (0,T),\\ & u^n(0)=u_0,\mathrm{in}\Omega .\end{array}$$

$$u^n\in C(\overline{\Omega }\times [0,T]),u_1\le u^n(x,t)\le u_3,(x,t)\in \overline{\Omega }\times [0,T].$$

$$u^n\in L^2(0,T;H^1(\Omega ))\cap C([0,T];L^2(\Omega )),u_t^n\in L^2(0,T;(H^1(\Omega ){)}^{\prime })+L^{4/3}(\Omega \times (0,T)),$$

$$\Vert u^n{\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2\le C,\Vert u^n{\Vert }_{L^2(0,T;H^1(\Omega ))}^2\le C,$$

$$\Vert u_t^n{|}_{L^{4/3}(0,T;(H^1(\Omega ){)}^{\prime })}^2\le C,\Vert {\varphi }_n\tilde{f}(u^n){\Vert }_{L^{4/3}(\Omega \times (0,T))}^2\le C,$$

$$\Vert u^{\epsilon }{\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2\le C,$$

$$u_1\le u^{\epsilon }(x,t)\le u_3,\mathrm{in}\overline{\Omega }\times [0,T].$$

$$\Vert w{\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}\le C|{\omega }_{\epsilon }{|}^{1/2},$$

$$\Vert w{\Vert }_{L^2(0,T;H^1(\Omega ))}\le C|{\omega }_{\epsilon }{|}^{1/2}.$$

$$\Vert w{\Vert }_{L^2(\Omega \times (0,T))}\le C|{\omega }_{\epsilon }{|}^{\frac{1}{2}+\beta }.$$

Full text

1 Introduction

Mathematical and numerical models of computational electrophysiology can provide quantitative tools to describe electrical heart function and disfunction [37], often complementing imaging techniques (such as computed tomography and magnetic resonance) for diagnostic and therapeutic purposes. In this context, detecting pathological conditions or reconstructing model features such as tissue conductivities from potential measurements yield to the solution of an inverse boundary value problem. Standard electrocardiographic techniques attempt to infer electrophysiological processes in the heart from body surface measurements of the electrical potential, as in the case of electrocardiograms (ECGs), or body surface ECGs (also known as body potential maps). These measurements can provide useful insights for the reconstruction of the cardiac electrical activity within the so-called electrocardiographic imaging, by solving the well-known inverse problem of electrocardiography111The inverse problem of electrocardiography aims at recovering the epicardial potential (that is, at the heart surface) from body surface measurements [36, 19, 18]. Since the torso is considered as a passive conductor, such an inverse problem involves the linear steady diffusion model as direct problem. A step further, aiming at computing the potential inside the heart from the epicardial potential, has been considered, e.g., in [11].. A much more invasive option to acquire potential measurements is represented by non-contact electrodes inside a heart cavity to record endocardial potentials.

Here we focus on the problem of detecting the position and the size of myocardial ischemias from a single boundary measurement of the electric potential. Ischemia is a reversible precursor of heart infarction caused by partial occlusion of one or more coronary arteries, which supply blood to the heart. If this condition persists, myocardial cells die and the ischemia eventually degenerates in infarction. For the time being, we consider an insulated heart model, neglecting the coupling with the torso; this results in the inverse problem of detecting inhomogeneities for a nonlinear parabolic reaction-diffusion equation (in our case, the so-called monodomain equation) dealing with a single measurement of the endocardial potential. Our long-term goal is indeed to deal with an inverse problem for the coupled heart-torso model, in order to detect ischemias from body surface measurements, such as those acquired on each patient with symptoms of cardiac disease through an ECG. The problem we consider in this paper is a mathematical challenge itself, almost never considered before. Indeed, difficulties include the nonlinearity of both the direct and the inverse problem, as well as the lack of measurements at disposal. Indeed, even for the linear counterpart of the inverse problem, it has been shown in [22] and [27] that infinitely many measurements are needed to detect uniquely the unknown inclusions, and that the continuous dependence of the inclusion from the data is logarithmic [20]. Moreover, despite the inverse problem of ischemia identification from measurements of surface potentials has been tackled in an optimization framework for numerical purposes [32, 30, 1, 15], a detailed mathematical analysis of this problem has never been performed. To our knowledge, no theoretical investigation of inverse problems related with ischemia detection involving the monodomain and/or the bidomain model has been carried out. On the other hand, recent results regarding both the analysis and the numerical approximation of this inverse problem in a much simpler stationary case have been obtained in [7, 8]. In order to obtain rigorous theoretical results additional assumptions are needed, for instance by considering small-size conductivity inhomogeneities. We thus model ischemic regions as small inclusions $\omega_{\epsilon}$ where the electric conductivity is significantly smaller than the one of healthy tissue and there is no ion transport.

We establish a rigorous asymptotic expansion of the boundary potential perturbation due to the presence of the inclusion adapting to the parabolic nonlinear case the approach introduced by Capdeboscq and Vogelius in [12] for the case of the linear conductivity equation. The theory of detection of small conductivity inhomogeneities from boundary measurements via asymptotic techniques has been developed in the last three decades in the framework of Electric Impedence Tomography (see, e.g., [5, 24, 13]). A similar approach has also been used in Thermal Imaging (see, e.g., [4]). We use these results to set a reconstruction procedure for detecting the inclusion. To this aim, as in [8], we propose a reconstruction algorithm based on topological optimization, where a suitable quadratic functional is minimized to detect the position and the size of the inclusion (see also [14]).

Numerical results obtained on an idealized left ventricle geometry assess the feasibility of the proposed procedure. Several numerical test cases also show the robustness of the reconstruction procedure with respect to measurement noise, unavoidable when dealing with real data. The modeling assumption on the small size of the inclusion, instrumental to the derivation of our theoretical results, is verifed in practice in the case of residual ischemias after myocardial infarction. On the other hand, a fundamental task of ECG’s imaging is to detect the presence of ischemias as precursor of heart infarction without any constraint on its size. For this reason, we also consider the case of the detection of larger size inclusions, for which the proposed algorithm still provides useful insights.

The paper is organized as follows. In Section 2 we describe the monodomain model of cardiac electrophysiology we are going to consider. In Section 3 we show some suitable wellposedness results concerning the direct problems, in the unperturbed (background) and perturbed cases. In Section 4 we prove useful energy estimates of the difference of the solutions of the two previous problems. The asymptotic expansion formula is derived in Section 5 and the reconstruction algorithm in Section 6. Numerical results are finally provided in Section 7. The appendix, Section 8, is devoted to a technical proof of a result needed in section 6.

2 The monodomain model of cardiac electrophysiology

The monodomain equation is a nonlinear parabolic reaction-diffusion PDE for the transmembrane potential, providing a mathematical description of the macroscopic electric activity of the heart [41, 18]. Throughout the paper we consider the following (background) initial and boundary value problem

[TABLE]

where $\Omega\subset{\bf R}^{3}$ is a bounded set with boundary $\partial\Omega$, and $k_{0}\in\mathbb{R},k_{0}>0$. Here $\Omega$ is the domain occupied by the ventricle, $u$ is the (transmembrane) electric potential, $f(u)$ is a nonlinear term modeling the ionic current flows across the membrane of cardiac cells, $k_{0}$ is the conductivity tensor of the healthy tissue, $C_{m}>0$ and $\nu>0$ are two constant coefficients representing the membrane capacitance and the surface area-to-volume ratio, respectively. For the sake of simplicity we deal with an insulated heart, namely we do not consider the effect of the surrounding torso, which behaves as a passive conductor. The initial datum $u_{0}$ represents the initial activation of the tissue, arising from the propagation of the electrical impulse in the cardiac conduction system. This equation yields a macroscopic model of the cardiac tissue, arising from the superposition of intra and extra cellular media, both assumed to occupy the whole heart volume (bidomain model), making the hypothesis that the extracellular and the intracellular conductivities are proportional quantities. Concerning the mathematical analysis of both the monodomain and the bidomain models, some results on the related direct problems have been obtained for instance in [6, 9, 10, 18].

We thus assume a phenomenological model to describe the effect of ionic currents through a nonlinear function of the potential. We neglect the coupling with the ODE system modeling the evolution of the so-called gating variables, which represent the amount of open channels per unit area of the cellular membrane and thus regulate the transmembrane currents.

In the case of a single gating variable $w$, a well-known option would be to replace $f$ by $g=g(u,w)$ where

[TABLE]

and $w$ solves the following ODE initial value problem, $\forall\,x\in\Omega$,

[TABLE]

for suitable (constant) parameters $\beta$, $\alpha$, $\rho$, $\gamma$. This is the so-called FitzHugh-Nagumo model for the ionic current, and the gating variable $w$ is indeed a recovery function allowing to take into account the depolarization phase. See, e.g., [18] for more details. In our case, the model (2.1) is indeed widely used to characterize the large-scale propagation of the front-like solution in the cardiac excitable medium.

As suggested in [18, Sect. 4.2] and [41, Sect. 2.2], hereon we consider the cubic function

[TABLE]

where $A>0$ is a parameter determining the rate of change of $u$ in the depolarization phase, and $u_{1}<u_{2}<u_{3}$ are given constant values representing the resting, threshold and peak potentials, respectively. Possible values of the parameters are, e.g., $u_{1}=-85mV$, $u_{2}=-65mV$ and $u_{3}=40mV$, $A=0.04$, see [41]. Note that both the sharpness of the wavefront and its propagation speed strongly depend on the value of the parameter $A$.

Consider now a small inhomogeneity located in a measurable bounded domain $\omega_{\varepsilon}\subset\Omega$, such that there exist a compact set $K_{0}$, with $\omega_{\varepsilon}\subset K_{0}\subset\Omega$, and a constant $d_{0}>0$ satisfying

[TABLE]

Moreover, we assume

[TABLE]

In the inhomogeneity $\omega_{\varepsilon}$ the conductivity coefficient and the nonlinearity take different values with respect the ones in $\Omega\backslash\omega_{\varepsilon}$. The problem we consider is therefore

[TABLE]

where $\chi_{D}$ stands for the characteristic function of a set $D\subset\mathbb{R}^{3}$. Here

[TABLE]

with $k_{0},k_{1}\in\mathbb{R},\,k_{0}>k_{1}>0$.

3 Well posedness of the direct problem

Problem (2.1) thus describes the propagation of the initial activation $u_{0}$ in an insulated heart portion (e.g., the left ventricle), and hereon will be referred to as the background problem; we devote Section 3.1 to the analysis of its well-posedness. The well-posedness of the perturbed problem modeling the presence of a small ischemic region in the domain will be instead analyzed in Section 3.2.

3.1 Well posedness of the background problem

For the sake of simplicity, throughout the paper we set $\nu=C_{m}=1$ and we assume that

[TABLE]

[TABLE]

Moreover, let us set

[TABLE]

The following well posedness result holds.

Theorem 3.1**.**

Assume (2.2), (3.1), (3.2). Then problem (2.1) admits a unique solution $u\in C^{2+\alpha,1+\alpha/2}(\overline{\Omega}\times[0,T])$ such that

[TABLE]

[TABLE]

where $C$ is a positive constant depending (at most) on $k_{0},T,\Omega,M_{1},M_{2},\|u_{0}\|_{C^{2+\alpha}(\overline{\Omega})}$.

*Proof. * We omit the details of the proof since (3.4) can be easily obtained using the results in [34, def. 3.1 and Thm. 4.1] and (3.5) by means of [29, Thm. 5.1.17 (ii) and Thm. 5.1.20].

$\Box$

3.2 Well posedness of the perturbed problem

The well-posedness of the perturbed problem (2.5) is provided by the following theorem.

Theorem 3.2**.**

Assume (2.2), (2.6), (3.1), (3.2). Then problem (2.5) admits a unique solution $u^{\varepsilon}$ such that

[TABLE]

Moreover, $u^{\varepsilon}\in C^{\alpha,\alpha/2}(\overline{\Omega}\times[0,T])$ and the following estimate holds

[TABLE]

where $C$ is a positive constant depending (at most) on $k_{0},k_{1},T,\Omega,\|u_{0}\|_{C^{\alpha}(\overline{\Omega})}$ and $M_{1}$.

*Proof. * Throughout the proof $C$ will be as in the statement of the Theorem.

Recalling the definition of $f$, there exist $k\geq 0,\,\alpha_{1}>0,\,\alpha_{2}>0,\,\lambda>0$ such that

[TABLE]

We formulate problem (2.5) in the weak form

[TABLE]

Setting $\tilde{f}(u)=f(u)-u$, (3.8) becomes

[TABLE]

Observe that, thanks to following the Poincaré type inequality in [7, formula (A.4)]

[TABLE]

the bilinear form $a_{\varepsilon}(u^{\varepsilon},v)=\left(\int_{\Omega}k_{\varepsilon}\nabla u^{\varepsilon}\cdot\nabla v+\int_{\Omega\backslash\omega_{\varepsilon}}u^{\varepsilon}v\right)$ is coercive. Indeed

[TABLE]

where $S$ is a positive constant depending on $\Omega$ and $k_{1}$.

Through the classical Faedo-Galerkin approximation scheme it is possible to prove that problem (2.5) admits a unique weak solution $u^{\varepsilon}$ satisfying (3.6).

In order to obtain further regularity for $u^{\varepsilon}$, let $\{\phi_{n}\}$ be a sequence such that

[TABLE]

and formulate the approximating problems

[TABLE]

Using the same arguments as in the proof of Theorem 3.1, we can prove that, $\forall\,n\in\mathbb{N}$, problem (3.12) admits a unique solution $u^{n}$ such that

[TABLE]

Moreover, by means again of a standard Faedo-Galerkin approximation scheme (for any $n$) we can prove that the solution to problem (3.12) satisfies

[TABLE]

[TABLE]

[TABLE]

where $C$ are some positive constants independent of $n$.

An application of [38, Thm. 8.1] implies that, up to a subsequence, $u^{n}\to\zeta$ strongly in $L^{2}(\Omega\times(0,T))$, so that $u^{n}\to\zeta$ a.e. in $\Omega\times(0,T)$ and $\phi_{n}\widetilde{f}(u^{n})\rightharpoonup\chi_{\omega_{\varepsilon}}\widetilde{f}(\zeta)\,\,{\rm in}\,\,L^{4/3}(\Omega\times(0,T))$. Since problem (2.5) has a unique solution (cf. (3.8)), we conclude that $\zeta=u^{\varepsilon}$ and satisfies

[TABLE]

Considering now the interior regularity result in [21, Theorem 2.1] (see also [28]) and the regularity up to the boundary contained in [21, Theorem 4.1], then we deduce (3.7). $\Box$

4 Energy estimates for $\mathbf{u^{\varepsilon}-u}$

In this section we prove some energy estimates for the difference between $u^{\varepsilon}$ and $u$, solutions to problem (2.5) and problem (2.1), respectively, that are crucial to establish the asymptotic formula for $u^{\varepsilon}-u$ of Theorem 4 in Section 5.

Proposition 4.1**.**

Assume (2.2), (2.6), (3.1), (3.2). Setting $w:=u^{\varepsilon}-u$, then

[TABLE]

[TABLE]

Moreover, there exists $0<\beta<1$ such that

[TABLE]

Here $C$ stands for a positive constant depending (at most) on $k_{0},k_{1},\Omega,T,M_{1},M_{2},\|u_{0}\|_{C^{2+\alpha}(\overline{\Omega})}$.

*Proof. * Throughout the proof $C$ will be as in the statement of the Theorem.

On account of the assumptions, Theorems 3.1 and 3.2 hold. Then $w$ solves the problem

[TABLE]

where we have set $\widetilde{k}:=k_{0}-k_{1}>0$ and

[TABLE]

$z_{\varepsilon}$ being a value between $u^{\varepsilon}$ and $u$. By means of (3.4), (3.13) and recalling (3.3), we have

[TABLE]

Multiplying the first equation by $w$ in (4.4) by $w$ and integrating by parts over $\Omega$, we get

[TABLE]

Adding and subtracting $\displaystyle{\int_{\Omega}\chi_{\Omega\setminus\omega_{\varepsilon}}w^{2}}$ and applying (3.11) we obtain

[TABLE]

Recalling (3.3) and (4.6), thanks to Young’s inequality we deduce

[TABLE]

so that

[TABLE]

and finally, see (3.5),

[TABLE]

Recalling $w(0)=0$, an application of Gronwall’s Lemma implies

[TABLE]

so that (4.1) follows. Integrating now inequality (4.7) on $(0,T)$ we get

[TABLE]

and a combination with (4.8) gives (4.2).

In order to obtain the more refined estimate (4.3), observe that $w$ also solves problem

[TABLE]

Let’s now introduce the auxiliary function $\overline{w}$, solution to the adjoint problem

[TABLE]

By the change of variable $t\to T-t$, problem (4.10) is equivalent to

[TABLE]

where we have set $z(x,t)=\overline{w}(x,T-t),\,\hat{p}_{\varepsilon}(x,t)=p_{\varepsilon}(x,T-t),\,\hat{w}(x,t)=w(x,T-t)$.

Since $|\chi_{\Omega/\omega_{\varepsilon}}\hat{p}_{\varepsilon}|$ is bounded in $\Omega\times(0,T)$ and $w\in C^{\alpha,\alpha/2}(\overline{\Omega}\times[0,T])$, standard arguments show that problem (4.11) admits a unique solution $z$ such that (see [28, Ch.4])

[TABLE]

Hereon, up to equation (4.15), all the equations depend on $t$ and are valid for every $t\in(0,T)$; however, we will omit this dependence for the sake of notation.

Moreover, multiplying the first equation in (4.11) by $z$ and integrating over $\Omega$, we get

[TABLE]

By means of Young’s inequality and recalling (4.6), we have

[TABLE]

and then

[TABLE]

Recalling $z(x,0)=0$, an application of Gronwall’s Lemma gives

[TABLE]

Let’s now multiply the first equation in (4.11) by $z_{t}$ and integrate over $\Omega$. We get

[TABLE]

An application of Young’s inequality gives

[TABLE]

and then

[TABLE]

Combining (4.15) and (4.14), integrating in time on $(0,t)$, and using $\nabla z(0)=0$ we deduce

[TABLE]

so that

[TABLE]

The same computations also gives

[TABLE]

Then, an application of standard elliptic regularity results to problem (4.11) implies (see [26])

[TABLE]

Recalling the definition of $z$ and $\hat{w}$, by estimates (4.16) and (4.18) we get

[TABLE]

Finally, we want to prove that there exists $p>2$ such that

[TABLE]

To this aim, on account of (4.19) and Sobolev immersion theorems, we deduce

[TABLE]

Moreover, again from (4.19) we have

[TABLE]

From well-known interpolation estimates (cf. [33]) we infer

[TABLE]

and therefore, using (4.19),

[TABLE]

so that (4.20) holds for any $p\in(2,\frac{10}{3}]$.

Let us now multiply the evolution equation in (4.9) by $\overline{w}$ and the evolution equation in (4.10) by $w$, respectively. Integrating on $\Omega$ we obtain

[TABLE]

[TABLE]

Summing up (4.25) and (4.26) we obtain

[TABLE]

subsequently, an integration in time on $(0,T)$ gives

[TABLE]

Recalling the conditions at time $t=0$ for $w$ and at time $t=T$ for $\overline{w}$, we get

[TABLE]

So that (4.27) becomes

[TABLE]

Using now Hölder inequality we deduce

[TABLE]

where we may choose for instance $p=10/3$ and $q=10/7$.

By means of (4.20) and (3.3), from the previous inequality we get

[TABLE]

and therefore

[TABLE]

Thanks to (3.5) we also have

[TABLE]

Finally, using again Hölder inequality and (4.2), and recalling that $q\in[10/7,)2$, we obtain

[TABLE]

Combining the previous estimate with (4.29), since $\frac{1}{q}\in(\frac{1}{2},\frac{7}{10}]$ we can conclude that (4.3) holds with $\beta\in(0,\frac{1}{5}]$. $\Box$

5 The asymptotic formula

In this section we derive and prove an asymptotic representation formula for $w=u_{\varepsilon}-u$ in analogy with [7] and [12]. Let $\Phi=\Phi(x,t)$ be any solution of

[TABLE]

Our main result is the following

Theorem 5.1**.**

Assume (2.2), (2.6), (3.1), (3.2). Let $u^{\varepsilon}$ and $u$ be the solutions to (2.5) and (2.1) and $\Phi$ a solution to (5.1), respectively. Then, there exist a sequence $\omega_{\varepsilon_{n}}$ satisfying (2.3) and (2.4) with $|\omega_{\varepsilon_{n}}|\rightarrow 0$, a regular Borel measure $\mu$ and a symmetric matrix $\mathcal{M}$ with elements $\mathcal{M}_{ij}\in L^{2}(\Omega,d\mu)$ such that, for $\varepsilon\to 0$,

[TABLE]

To prove Theorem 5.1, we need to state some preliminary results. Let $v_{\varepsilon}^{(j)}$ and $v^{(j)}$ be the variational solutions (depending only on $x\in\Omega$) to the problems

[TABLE]

$n_{j}$ being the $j-th$ coordinate of the outward normal to $\partial\Omega$. It can be easily verified that

[TABLE]

The following results hold

Proposition 5.2**.**

Let $v_{\varepsilon}^{(j)}$ and $v^{(j)}$ solutions to (5.3), then there exists $C(\Omega)>0$ such that

[TABLE]

Moreover, for some $\eta\in(0,\frac{1}{2})$, there exists $C(\Omega,\eta)>0$ such that

[TABLE]

*Proof. * See Lemma 1 in [12]. $\Box$

Proposition 5.3**.**

Let $u$ and $u_{\varepsilon}$ be the solutions to problems (2.1) and $\eqref{perturbed}$, respectively. Consider $v_{\varepsilon}^{(j)}$ and $v^{(j)}$ as in (5.3). Then, for any $\Phi\in C^{1}(\overline{\Omega}\times[0,T])$ with $\Phi(x,T)=0$, the folllowing holds as $\varepsilon\to 0$,

[TABLE]

*Proof. * We follow the ideas in [7] and [12]. Since $w=u_{\varepsilon}-u$, then we obtain the identity

[TABLE]

Moreover, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A combination with (5.8) gives

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Then, on account of (4.1), (4.2), (5.5), (5.6) and Schwarz inequality, we get

[TABLE]

Let us consider now problem (4.4). Multiplying both sides of the first equation by $v_{\varepsilon}^{(j)}\Phi$, and integrating by parts on $\Omega\times(0,T)$, we obtain

[TABLE]

On the other hand, multiplying the first equation in (4.9) by $v^{(j)}\Phi$ and integrating by parts on $\Omega\times(0,T)$, we get

[TABLE]

A combination of (5.9), (5.10) and (5.11) gives, for $\varepsilon\to 0$,

[TABLE]

[TABLE]

[TABLE]

from which we deduce

[TABLE]

[TABLE]

By means of (4.29), (4.3), (5.5) and (5.6), and recalling also (3.3) and (4.6), an application of the Hölder inequality both in space and time gives, for $\varepsilon\to 0$,

[TABLE]

and then

[TABLE]

Consider the first term in the last line of (5.12). Integrating by parts in time and recalling that $\Phi(T)=0$, $w(0)=0$, $(v^{(j)}-v_{\varepsilon}^{(j)})_{t}=0$, we finally have (cf. also (4.3), (5.6)), for $\varepsilon\to 0$,

[TABLE]

Combining (5.12) and (5.13) we get

[TABLE]

then formula (5.7) is true. $\Box$

Proof of Theorem 5.1. Following [12, Sec.3], there exist a regular Borel measure $\mu$, a symmetric matrix $\mathcal{M}$ with elements $\mathcal{M}_{ij}\in L^{2}(\Omega,d\mu)$, a sequence $\omega_{\varepsilon_{n}}$ with $|\omega_{\varepsilon_{n}}|\to 0$ such that

[TABLE]

in the weak∗ topology of $C^{0}(\overline{\Omega})$. On account of (3.5), we deduce also

[TABLE]

in the weak∗ topology of $C^{0}(\overline{\Omega})$. Moreover, recalling (3.5), (4.3) and (5.4), we get

[TABLE]

where $C$ is independent of $\varepsilon_{n}$. Hence

[TABLE]

in the weak∗ topology of $C^{0}(\overline{\Omega}\times[0,T])$. Combining (5.7), (5.16) and (5.18) we obtain

[TABLE]

Now multiply the first equation in (5.1) by $w$ and the first equation in (4.9) by $\Phi$ on $\Omega\times(0,T)$.

Then, integrating by parts, we get

[TABLE]

[TABLE]

Summing up the two previous equations, we have

[TABLE]

Observe that the following identities hold

[TABLE]

and then, from (5.20) we infer

[TABLE]

Moreover, on account of (4.3), we have

[TABLE]

The last equality in (5.23) is a consequence of the regularity of $f$ (see (4.5), from which $|p_{\varepsilon}-f^{\prime}(u)|\leq C|w|$ follows) and (4.3). Combining (5.22) and (5.23) we obtain

[TABLE]

And finally, by means of (5.15), (5.18) and (5.19), the formula (5.2) holds. $\Box$

Remark 5.4**.**

We would like to emphasize that, with minor changes, the asymptotic expansion extends to the case of piecewise smooth anisotropic conductivities of the form

[TABLE]

where $\mathbb{K}_{0},\mathbb{K}_{1}\in C^{\infty}(\Omega)$ are symmetric matrix valued functions satisfying

[TABLE]

with $0<\alpha_{1}<\beta_{1}<\alpha_{0}<\beta_{0}$. Then, the asymptotic formula (5.2) becomes

[TABLE]

where $\Phi$ solves

[TABLE]

and $u$ is the background solution of

[TABLE]

The matrix $\mathcal{M}$ is called the polarization tensor associated to the inhomogeneity $\omega_{\varepsilon}$. Indeed, all the results of the previous sections can be extended to the case of constant anisotropic coefficients using for instance the regularity results contained in [28].

6 A reconstruction algorithm

We now use the asymptotic expansion derived in the previous section to set a reconstruction procedure for the inverse problem of detecting a spherical inhomogeneity $\omega_{\varepsilon}$ from boundary measurements of the potential. Following the approach of [8], [14], but taking into account the time-dependence of the problem, we introduce the mismatch functional

[TABLE]

being $u^{\varepsilon}$ the solution of the perturbed problem (2.5) in presence of an inclusion $\omega_{\varepsilon}$ satisfying hypotheses (2.3), (2.4). It is possible to reformulate the inverse problem in terms of the following minimization problem

[TABLE]

among all the small inclusions, well separated from the boundary. We introduce the following additional assumption on the exact inclusion

[TABLE]

being $z\in\Omega$ and $B$ an open, bounded, regular set containing the origin. We remark that we prescribe the geometry of the inclusion to be fixed throughout the whole observation time. The restriction of the functional $J$ to the class of inclusions satisfying (6.3) is denoted by $j(\varepsilon;z)$. We can now define the topological gradient $G:\Omega\rightarrow\mathbb{R}$ as the first order term appearing in the asymptotic expansion of the cost functional with respect to $\varepsilon$, namely

[TABLE]

where $j(0)=\int_{0}^{T}\int_{\partial\Omega}(u-u_{meas})^{2}$ and $u$ is the solution of the unperturbed problem (2.1). Under the assumption that the exact inclusion has small size and satisfies hypothesis (6.3), a reconstruction procedure consists in identifying the point $\bar{z}\in\Omega$ where the topological gradient $G$ attaints its minimum. Indeed, the cost functional achieves the smallest value when it is evaluated in the center of the exact inclusion. Thanks to the hypothesis of small size, we expect the reduction of the cost functional $j$ to be correctly described by the first order term $G$, up to a reminder which is negligible with respect to $\varepsilon$.

Nevertheless, in order to define a reconstruction algorithm, we need to efficiently evaluate the topological gradient $G$. According to the definition,

[TABLE]

Evaluating $G$ in a single point $z\in\Omega$ would require to solve the direct problem several times in presence of inclusions centered at $z$ with decreasing volume. This procedure can be indeed avoided thanks to a useful representation formula that can be deduced from the asymptotic expansion (5.2). To show this we need the following preliminary Proposition the proof of which is given in the Appendix.

Proposition 6.1**.**

Consider the problem

[TABLE]

Given a compact set $K\subset\Omega$ such that $d(K,\partial\Omega)\geq d_{0}>0$ the following estimate holds

[TABLE]

On account of Proposition 6.1, we deduce the following representation of the topological gradient

Proposition 6.2**.**

The topological gradient of the cost functional $j(\varepsilon,z)$ can be expressed by

[TABLE]

where $W$ is the solution of the adjoint problem:

[TABLE]

*Proof. * Consider the difference

[TABLE]

According to (5.2) and to the definition of the adjoint problem (6.7), we can express

[TABLE]

Since we assume (6.3), the measure $\mu$ associated to the inclusion is the Dirac mass $\delta_{z}$ centered in point $z$ (see [12]). Hence

[TABLE]

Moreover, by (5.2), the second term in the left-hand side of (6.8) can be expressed as

[TABLE]

where $\Phi$ is the solution to (6.4). Thanks to regularity results on $u$ (see Theorem 3.1) and using Proposition 6.1 with $K=\Omega_{d_{0}}=\{x\in\Omega\textit{ s.t. }d(x,\partial\Omega)\geq d_{0}\}$, we obtain

[TABLE]

Replacing (6.9) and (6.10) in (6.8), we finally get

[TABLE]

$\Box$

Thanks to the representation formula (6.6), evaluating the topological gradient of the cost functional requires just the solution of two initial and boundary value problems. This yields the definition of a one-shot algorithm for the identification of the center of a small inclusion satisfying hypotesis (6.3) (see Algorithm 1).

Inspired by the electrophysiological application, we consider moreover the possibility to have partial boundary measurements, i.e. the case where the support of the function $u_{meas}$ is not the whole boundary $\partial\Omega$ but only a subset $\Gamma\subset\partial\Omega$. In this case, it is possible to formulate a slightly different optimization problem, in which we aim at minimizing the mismatch between the measured and the perturbed data just on the portion $\Gamma$ of the boundary. The same reconstruction algorithm can be devised for this problem, by simply changing the expression of the Neumann condition of the adjoint problem.

7 Numerical results

In order to implement Algorithm 1 for the detection of inclusions, it is necessary to approximate the solution of the background problem (2.1) and the adjoint problem (6.7). Moreover, when considering synthetic data $u_{meas}$, we must be able to compute the solution to the perturbed problem (2.5) in presence of the exact inclusion. We rely on the Galerkin finite element method for the numerical approximation of these problems. The one-shot procedure makes the reconstruction algorithm very efficient, only requiring the solution of an adjoint problem for each acquired measurement over the time interval, without entailing any iterative (e.g. descent) method for numerical optimization.

7.1 Finite Element approximation of initial and boundary value problems

The background problem (2.1) can be cast in weak form as follows

$\forall\ t\in(0,T)$, find $u(t)\in V=H^{1}(\Omega)$ such that $u(0)=u_{0}$ and

[TABLE]

By introducing a finite-dimensional subspace $V_{h}$ of $V$, $dim(V_{h})=N_{h}<\infty$, the Galerkin (semi-discretized in space) formulation of problem (7.1) reads

$\forall\,t\in(0,T)$, find $u_{h}(t)\in V_{h}$ such that $u_{h}(0)=u_{h,0}$ and

[TABLE]

where $(\cdot,\cdot)$ is the inner product in $L^{2}(\Omega)$, $b(u,v)=\int_{\Omega}k_{1}\nabla u\cdot\nabla v$, $F(u,v)=\int_{\Omega}f(u)v$, $f$ is defined as in (2.2) and $u_{h,0}$ is the $H^{1}$-projection of $u_{0}$ onto $V_{h}$.

To obtain a full discretization of the problem, we introduce a finite difference approximation in time. According to the strategy reported in [18], [16], we rely on a semi-implicit scheme which allows an efficient treatment of the nonlinear terms. Let us consider an uniform partition $\{t^{n}\}_{n=0}^{N}$ of the time interval $[0,T]$ of step $\tau=\frac{T}{N}$ s.t. $t^{0}=0,\ t^{N}=T$. Then, the fully discrete formulation of (2.1) is given by

$\forall\,n=0,\ldots N-1$, find $u_{h}^{n+1}\in V_{h}$ such that $u_{h}^{0}=u_{0,h}$ and

[TABLE]

With the same discretization strategy one may describe a numerical scheme for the approximate solution of the perturbed problem, using the weak form reported in (3.8) and introducing the forms

[TABLE]

The adjoint problem, instead, requires the introduction of the form $dF(u,v;w)=\int_{\Omega}f^{\prime}(w)uv$, which is bilinear with respect to $u$ and $v$. Thanks to the linearity of the adjoint problem, we can consider a fully implicit Crank-Nicolson scheme

$\forall\,n=0,\ldots N-1$, find $w_{h}^{n}\in V=H^{1}(\Omega)$ such that $w^{N}_{h}=0$ and

[TABLE]

The existence and the uniqueness of the solutions of the fully-discrete problems (7.3) and (7.4) follow by the well-posedness of the continuous problems, since $V_{h}$ is a subspace of $H^{1}(\Omega)$. For further details on the stability and of the convergence of the proposed schemes we refer to [23], [40] and [18].

The numerical setup for the simulation is represented in Figure 1. We consider an idealized geometry of the left ventricle (which has been object of several studies, see e.g. [18], [16]), and define a tetrahedral tesselation $\mathcal{T}_{h}$ of the domain. The discrete space $V_{h}$ is the P1-Finite Element space over $\mathcal{T}_{h}$, i.e. the space of the continuous functions over $\Omega$ which are linear polynomials when restricted on each element $T\in\mathcal{T}_{h}$. The mesh we use for all the reported results consists of 24924 tetrahedric elements and $N_{h}=5639$ nodes. We report also the anisotropic structure considered in all the recontruction tests, according to [39] and [18]. The conductivity matrix $\mathbb{K}_{0}$ for the monodomain equation is given by $\mathbb{K}_{0}(x)=\mathbb{K}^{e}(x)(\mathbb{K}^{e}(x)+\mathbb{K}^{i}(x))^{-1}\mathbb{K}^{i}(x)$, where both $\mathbb{K}^{i}$ and $\mathbb{K}^{e}$ are orthotropic tensors with three constant positive real eigenvalues, namely

[TABLE]

The eigenvectors $\vec{e_{f}}$, $\vec{e_{t}}$ and $\vec{e_{r}}$ are associated to the three principal directions of conductivity in the heart tissue: respectively, the fiber centerline, the tangent direction to the heart sheets and the transmural direction (normal to the sheets).

For the direct problem simulations, we consider the formulation reported in (2.1), specifying realistic values for the parameters $C_{m}$ and $\nu$. We have rescaled the values of the coefficients $u_{1}$, $u_{2}$, $u_{3}$ and $A^{2}$ in order to simulate the electric potential in the adimensional range $[0,1]$. The rescaling is performed by the transformation $\tilde{u}=(\alpha+u)/\beta$, where $\alpha=0.085mV$ and $\beta=0.125mV$, whereas for the sake of simplicity we will still denote by u the rescaled variable $\tilde{u}$. We consider the initial datum $u_{0}$ to be positive on a band of the endocardium, representing the initial stimulus provided by the heart conducting system. The most important parameters, considered in accordance with [25], [41], are reported in Table 1.

In Figure 2 we report the solution of the discrete background problem (7.3) at different time instants, comparing the isotropic and the anisotropic cases.

7.2 Reconstruction of small inclusions

We now tackle the problem of reconstructing the position of a small inhomogeneity using the knowledge of the electric potential of the tissue on a portion $\Gamma$ of the boundary of the domain. In particular, we assume that $u_{meas}$ is known on the endocardium, i.e. the inner surface of the heart cavity. We generate synthetic data on a more refined mesh and test the effectiveness of Algorithm 1 in the reconstruction of a small spherical inclusion in different positions. In Figure 3 we report the value assumed by the topological gradient, and superimpose the exact inclusions: we observe a negative region in proximity of the position of the real inclusion. The algorithm precisely identifies the region where the inclusion is present, whereas the minimum may in general be found along the endocardium also when the center of the real inclusion is not located on the heart surface. Nevertheless, due to the thinness of the domain the reconstructed position is close to the real one.

This slight loss in accuracy seems to be an intrinsic limit of the topological gradient strategy applied to the considered problem. We point out that the reconstruction is performed by relying on a single measurement acquired on the boundary, which is indeed a constraint imposed by the physical problem. Hence, all the techniques relying on the introduction of several measurements to increase the quality of the reconstruction are impracticable. A different strategy, as proposed in several works addressing the steady-state case, may consist in introducing a modification to the cost functional $J$. In [3] and many related works the authors introduce a cost functional inherited from imaging techniques, whereas in [14], [35] different strategies involving the Kohn-Vogelius functional or similar ones are explored. Nevertheless, the nonlinearity of the direct problem considered in this work prevents the possibility to apply these techniques, since the analitycal expressions of the fundamental solution, single and double layer potentials would not be available.

7.3 Reconstruction in presence of experimental noise

We test the stability of the algorithm in presence of experimental noise on the measured data $u_{meas}$. We consider different noise levels, according to the formula

[TABLE]

where $\eta(x,t)$, for each point $x$ and instant $t$, is a Gaussian random variable with zero mean and standard deviation equal to $u_{3}-u_{1}$, whereas $p\in[0,1]$ is the noise level. In Figure 4 the results of the reconstruction with different noise levels are compared. The algorithm shows to be highly stable with respect to high rates of noise, with increasing accuracy as the noise level reduces.

7.4 Reconstruction from partial discrete data

A further test case we have performed deals with the reconstruction of the position of small inclusions starting from the knowledge of partial data. We are interested in assessing the effectiveness of our algorithm when the electric potential is measured only on a discrete set of points on the endocardium, possibly simulating the procedure of intracavitary electric measurements. Figure 5 shows that the algorithm is able to detect the region where the small ischemia is located from the knowledge of the potential on $N_{p}=246,61,15$ different points. The position of the reconstructed inclusion is slightly affected by the reduction of sampling points; nevertheless, reliable reconstructions can be obtained even with a very small number of points.

For the same purpose, we have tested the capability of the reconstruction procedure to avoid false positives: the algorithm is able to distinguish the presence of a real ischemia from the case where no ischemic region is present, also in the case where the data are recovered only at a finite set of points, and are affected by noise. We compare the value of the cost functional and of the minimum of the topological gradient obtained through Algorithm 1 on data generated when $(i)$ a small ischemia is present in the tissue or $(ii)$ no inclusion is considered. The measurement is performed on a set of $N_{p}=100$ points and is affected by different noise levels. The results are reported in Table 2.

We observe that the presence of a small noise on the measured data causes a great increase of the cost functional $J$: with $5\%$ noise, e.g., the value of $J$ is two orders of magnitude greater than the value that $J$ assumes in presence of a small inclusion without noise. Nevertheless, the topological gradient $G$ allows to distinguish the false positive cases, since (at least in case of a small noise level) the value attained by its minimum in presence of a small inclusion is considerably lower than the random oscillations of $G$ due to the noise.

7.5 Reconstruction of larger inclusions

We finally assess the performance of Algorithm 1, developed for the reconstruction of small inclusions well separated from the boundary, in detecting the position of extended inclusions, which may be of greater interest in view of the problem of detecting ischemic regions. Indeed, total occlusion of a major coronary artery generally causes the entire thickness of the ventricular wall to become ischemic (transmural ischemia) or, alternatively, a significant ischemia only in the endocardium, that is, the inner layer of the myocardium (subendocardial ischemia). See, e.g., [17] for a detailed investigation of the interaction between the presence of moderate or severe subendocardial ischemic regions and the anisotropic structure of the cardiac muscle.

The most important assumption on which our one-shot procedure relies is that the variation of the cost functional from the value $J(0)$ attained in the background case to the value related to the exact inclusion can be correctly described by the first order term of its asymptotic expansion, the topological gradient $G$. Removing the hypothesis of the small size, we cannot rigorously assess the accuracy of the algorithm, however it still allows us to identify the location of the ischemic region.

The results reported in Figure 6 show that in presence of a inclusion of larger size (and not even separated from the boundary), the minimum of the topological gradient is found to be close to the position of the inclusion, and attains lower values with respect to the previously reported cases.

Moreover, in Figure 7 we also assess the stability of the reconstruction with respect to the presence of noisy data and partial measurements, as done in the case of small inclusions.

8 Conclusions and perspectives

A rigorous theoretical analysis of the inverse problem of detecting inhomogeneities in the monodomain equations has allowed us to set up a numerical reconstruction procedure, aiming at the detection of ischemic regions in the myocardic tissue from a single measurement of the endocardial potential. The identification is made possible by evaluating the topological gradient of a quadratic cost functional, requiring the solution of two initial and boundary value problems, the background problem and the adjoint one. Numerical results are encouraging and allow to estimate the position of the inclusion, although the identified inhomogeneity is nearly always detected on the boundary where the measurement is acquired. Nevertheless, provided a single measurement can be used for the sake of identification, and a one-shot procedure is performed, the obtained results give useful insights.

Many issues are still open. Concerning the mathematical model, an even more interesting case would be the one involving the heart-torso coupling is considered, so that more realistic (and noninvasive) body surface measurements can be employed. Setting and analyzing the inverse problem in this context represents the natural continuation of the present work. To close the gap between the rigorous mathematical setting and the practice, the two assumptions made in this work about the size of the inclusion and its distance from the boundary should be relaxed. Numerical results shown in Section 7.5 provide a first insight on the detection of inclusions with larger size, as those corresponding to transmural or subendocardial ischemias. From a mathematical standpoint, this problem is still open. Also in the case of a linear direct problem, very few results can be found in literature, see, e.g. [2]. Estimating the size of the inclusion is another open question in the case of parabolic PDEs, also for linear equations. The case of multiple inclusions, addressed in [8] for a stationary nonlinear problem, could also be considered. Last, but not least, the topological optimization framework addressed in this paper could also be combined with an iterative algorithm, such as the level set method, or with the solution of a successive shape optimization problem, to achieve a full reconstruction both of the dimension and the shape of the inclusion.

9 Appendix - Proof of Proposition 6.1

*Proof. * Setting $Z(t)=\Phi(T-t),\,t\in(0,T)$, we get an equivalent problem to (6.4)

[TABLE]

We will prove that $Z\in L^{2}(0,T;H^{3}(K))\hookrightarrow L^{1}(0,T;W^{1,\infty}(K))$. To this aim we need to derive some energy estimates. Multiplying the first equation in (9.1) by $Z$, an application of Young’s inequality leads to

[TABLE]

where $C=C(k_{0},M_{2},\Omega)>0$. An application of Gronwall’s lemma gives

[TABLE]

so that

[TABLE]

Instead, integrating (9.2) in time over $[0,t]$ we get

[TABLE]

and finally

[TABLE]

where $C$ is a positive constant depending on $k_{0},M_{2},\Omega,T$. We remark that, by standard regularity results, $Z$ is smooth on $E\times[0,T]$, for any compact $E\subset\Omega$.

Consider now two compact sets $K_{1}$ and $K_{2}$ such that

[TABLE]

It is possible to construct two functions $\xi_{1},\,\xi_{2}$ and two constants $b_{1},\,b_{2}$ satisfying

[TABLE]

[TABLE]

Let us multiply the first equation of (6.4) by $-\Delta Z$. Then the following holds

[TABLE]

Multiplying (9.5) by $\xi_{1}$, integrating on $\Omega\times(0,T)$ and using the definitions of $Z$ we get

[TABLE]

Combining (9.6) and the first equation in (9.1), applying Young’s inequality and taking into account (4.6) and the fact that $0\leq\xi\leq 1$, we obtain

[TABLE]

Integrating by parts the term $\int_{0}^{T}\int_{\Omega}\Delta Z\nabla Z\cdot\nabla\xi_{1}$, we can easily deduce

[TABLE]

where $C$ is a positive constant depending on $M_{2}$, $k_{0}$, $\xi_{1}$. Hence, since $\xi_{1}=1$ in $k_{0}$, we get

[TABLE]

Observe that, replacing $T$ by $t\in(0,T]$ in (9.7), we deduce also

[TABLE]

Combining (9.3), (9.8) and (9.9), we obtain

[TABLE]

where $C$ is a positive constant depending on $k_{0},M_{2},\Omega,T,\xi_{1}$.

On account of the first equation in (9.1) and the previous estimates, we get

[TABLE]

where $C$ is a positive constant depending on $k_{0},M_{2},\Omega,T,\xi_{1}$.

Now, let us multiply the first equation of (9.1) by $-\Delta Z_{t}$. We obtain

[TABLE]

Multiplying the previous equation by $\xi_{2}$ and integrating on $\Omega\times(0,T)$, then a suitable integration by parts in space implies

[TABLE]

[TABLE]

Integrating by parts the second term of the left-hand side and by parts in space the terms in the right-hand side, by an application of Young’s inequality we finally get

[TABLE]

where the constant $C>0$ depends on $\xi_{2},M_{2}$. A combination with (9.3), (9.4), (9.11) gives

[TABLE]

where the constant $C>0$ depends on $k_{0},M_{2},\Omega,T,\xi_{1},\xi_{2}$. In order to prove the desired regularity, we need to take into account also the third-order derivatives, in particular the operator $\nabla\Delta Z$. Observe that from the first equation in (9.1) we get

[TABLE]

Hence, we can conclude

[TABLE]

where $C$ is a positive constant depending on $k_{0},\frac{1}{k_{0}},M_{2},\Omega,T,\xi_{1},\xi_{2}$.

Recalling (9.10) and the fact that $K\subset K_{2}\subset K_{1}$, standard regularity results imply

[TABLE]

Finally, from (9.10) and (9.13), by Sobolev immersion theorems, we get

[TABLE]

where $C$ is a positive constant depending on $k_{0},\frac{1}{k_{0}},M_{2},\Omega,T,\xi_{1},\xi_{2}$.

Recalling the relation between $Z$ and $\Phi$ we get (6.5). $\Box$

Acknowledgments

E. Beretta, C. Cavaterra, M.C. Cerutti and L. Ratti thank the New York University in Abu Dhabi for its kind hospitality that permitted a further development of the present research. The work of C. Cavaterra was supported by the FP7-IDEAS-ERC-StG 256872 (EntroPhase) and by GNAMPA (Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni).