On the determination of ischemic regions in the monodomain model of cardiac electrophysiology from boundary measurements

TL;DR

This paper develops a mathematical method to detect small ischemic regions in the heart tissue by analyzing boundary measurements within the monodomain model of cardiac electrophysiology, using asymptotic expansions and topological gradients.

Contribution

It introduces a novel approach combining asymptotic analysis and topological derivatives to identify ischemic regions from boundary data in cardiac models.

Findings

The method accurately detects small ischemic regions.

Numerical experiments confirm the robustness of the detection algorithm.

The approach is theoretically grounded and computationally validated.

Abstract

In this paper we consider the monodomain model of cardiac electrophysiology. After an analysis of the well-posedness of the model, we determine an asymptotic expansion of the perturbed potential due to the presence of small conductivity inhomogeneities (modeling small ischemic regions in the cardiac tissue) and use it to detect the anomalies from partial boundary measurements. This is done by determining the topological gradient of a suitable boundary misfit functional. The robustness of the algorithm is confirmed by several numerical experiments.

Figures

Tables

Table 1. Table 1: Values of the main parameters

$k_{0,1}$	$k_{0,2}$	$k_{1,1}$	$k_{1,2}$	$A$	$a$	$ϵ$
$1.200$	$0.2538$	$0.2308$	$0.0062$	$8$	$0.15$	$0.05$

Equations

$$\{\begin{array}{rlr}{\partial }_{t}u-div(K_0\nabla u)+f(u,w)& =0& \text{in }{Q}_T,\\ K_0\nabla u\cdot \nu & =0& \text{on }{\Gamma }_T,\\ {\partial }_{t}w+g(u,w)& =0& \text{in }{Q}_T,\\ u(\cdot ,0)=u_{0}w(\cdot ,0)& =w_0& \text{in }\Omega ,\end{array}$$

$$\{\begin{array}{rlr}{\partial }_t{u}_{\omega }-div(K_{\omega }\nabla u_{\omega })+(1-{\chi }_{\omega })f(u_{\omega },w_{\omega })& =0& \text{in }{Q}_T,\\ K_{\omega }\nabla u_{\omega }\cdot \nu & =0& \text{on }{\Gamma }_T,\\ {\partial }_t{w}_{\omega }+g(u_{\omega },w_{\omega })& =0& \text{in }{Q}_T,\\ u_{\omega }(\cdot ,0)=u_0{w}_{\omega }(\cdot ,0)& =w_0& \text{in }\Omega ,\end{array}$$

$$\text{vspace}-0.25cm\exists {\mathcal{C}}_0\text{ compact subset of }\Omega \text{ s.t. }\omega \subset {\mathcal{C}}_0\text{ and }dist({\mathcal{C}}_0,\partial \Omega )\ge d_0>0.\text{vspace}-0.25cm$$

$$\text{vspace}-0.25cmf(u,w)=Au(u-a)(u-1)+uwg(u,w)=ϵ(Au(u-1-a)+w),$$

$$\vec{p}({\xi }_1,{\xi }_2)\cdot \left(\begin{array}{r}-f({\xi }_1,{\xi }_2)\\ -g({\xi }_1,{\xi }_2)\end{array}\right)\le 0\forall ({\xi }_1,{\xi }_2)\in \partial S.$$

$$⟨{\partial }_t{u}_{\omega },\phi {⟩}_{\ast }+{\int }_{\Omega }{K}_{\omega }\nabla u_{\omega }\cdot \nabla \phi +{\int }_{\Omega }(1-{\chi }_{\omega })f(u_{\omega },w_{\omega })\phi$$

$${\int }_{\Omega }{\partial }_t{w}_{\omega }\psi +{\int }_{\Omega }g(u_{\omega },w_{\omega })\psi$$

$$\exists \{{\varphi }_k\}\subset C_c^2(\Omega ):0\le {\varphi }_k(x)\le 1\forall x\in \Omega ,{\varphi }_k\to {\chi }_{\omega }\text{ in L2(Ω) and a.e.},$$

$$\{\begin{array}{rlr}{\partial }_t{u}_k-div(K_k\nabla u_k)+(1-{\varphi }_k)f(u_k,w_k)& =0& \text{in }{Q}_T,\\ K_0\nabla u_k\cdot \nu & =0& \text{on }{\Gamma }_T,\\ {\partial }_t{w}_k+g(u_k,w_k)& =0& \text{in }{Q}_T,\\ u_k(\cdot ,0)=u_0{w}_k(\cdot ,0)& =w_0& \text{in }\Omega ,\end{array}$$

$$\frac{1}{2}\frac{d}{dt}\left({\Vert u_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2+{\Vert w_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2\right)+{\int }_{\Omega }{K}_k\nabla u_k(\cdot ,t)\cdot \nabla u_k(\cdot ,t)$$

$$=-{\int }_{\Omega }(1-{\varphi }_k)f(u_k(\cdot ,t),w_k(\cdot ,t))u_k(\cdot ,t)-{\int }_{\Omega }g(u_k(\cdot ,t),w_k(\cdot ,t))w_k(\cdot ,t).$$

$$\frac{1}{2}\frac{d}{dt}\left({\Vert u_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2+{\Vert w_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2\right)+k_{min}{\Vert \nabla u_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2$$

$$\le \frac{1}{2}\left({\Vert u_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2+{\Vert w_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2\right)+\frac{1}{2}|\Omega |(M_f^2+M_g^2).$$

$$\left({\Vert u_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2+{\Vert w_k(\cdot ,t)\Vert }_{L^2(\Omega )}^2\right)$$

$$\le \left({\Vert u_0\Vert }_{L^2(\Omega )}^2+{\Vert w_0\Vert }_{L^2(\Omega )}^2+|\Omega |(M_f^2+M_g^2)t\right)e^t,$$

$${\Vert u_k\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2+{\Vert w_k\Vert }_{L^{\infty }(0,T;L^2(\Omega ))}^2$$

$$\le \left({\Vert u_0\Vert }_{L^2(\Omega )}^2+{\Vert w_0\Vert }_{L^2(\Omega )}^2+|\Omega |(M_f^2+M_g^2)T\right)e^T:=c_1^2.$$

$${\Vert u_k\Vert }_{L^2(0,T,H^1(\Omega ))}^2\le c_1^{2}T+\frac{1}{2k_{min}}\left((|\Omega |(M_f^2+M_g^2)+c_1^2)T+{\Vert u_0\Vert }_{L^2(\Omega )}^2+{\Vert w_0\Vert }_{L^2(\Omega )}^2\right)=:c_2^2.$$

$$\left|⟨{\partial }_t{u}_k(\cdot ,t),\phi {⟩}_{\ast }\right|$$

$${\Vert {\partial }_t{u}_k\Vert }_{L^2(0,T,H^{\ast })}\le (k_{max}{c}_2+M_f|\Omega {|}^{\frac{1}{2}})=:c_3.$$

$$u_k\stackrel{L^2(0,T,H^1)}{\rightharpoonup }u,{\partial }_t{u}_k\stackrel{L^2(0,T,H^{\ast })}{\rightharpoonup }{u}^{\ast },w_k\stackrel{L^2(Q_T)}{\rightharpoonup }w,{\partial }_t{w}_k\stackrel{L^2(Q_T)}{\rightharpoonup }{w}^{\ast }.$$

$$w_k(x,t)=e^{-ϵt}{w}_0(x)+ϵA{e}^{-ϵt}{\int }_0^t((1+a)u_k-u_k^2)e^{ϵs}ds.$$

$$⟨{\partial }_t{u}_k,\phi {⟩}_{\ast }+{\int }_{\Omega }{K}_k\nabla u_k\cdot \nabla \phi +{\int }_{\Omega }(1-{\varphi }_k)A{u}_k(u_k-a)(u_k-1)\phi$$

$$+{\int }_{\Omega }(1-{\varphi }_k)u_k{w}_k\phi +{\int }_{\Omega }{\partial }_t{w}_k\psi +{\int }_{\Omega }Aϵ(u_k^2-(1+a)u_k)\psi +{\int }_{\Omega }ϵw_k\psi =0.$$

$$⟨{\partial }_t{u}_{\omega },\phi {⟩}_{\ast }+{\int }_{\Omega }{K}_{\omega }\nabla u_{\omega }\cdot \nabla \phi +{\int }_{\Omega }(1-{\chi }_{\omega })f(u_{\omega },w_{\omega })\phi$$

$$+{\int }_{\Omega }{\partial }_t{w}_{\omega }\psi +{\int }_{\Omega }g(u_{\omega },w_{\omega })\psi =0.$$

$${\partial }_t{u}_{\omega }-div(K_{\omega }\nabla u_{\omega })=-(1-{\chi }_{\omega })f(u_{\omega },w_{\omega })$$

$$\text{find ω⊂Ω s.t. }{u}_{\omega }{|}_{\partial \Omega }=u_{meas}.$$

$$|{\omega }_{\epsilon }|\to 0\text{ as }\epsilon \to 0.$$

$${\int }_{\Omega }{\partial }_t{U}_{\epsilon }\phi +{\int }_{\Omega }{K}_{\epsilon }\nabla U_{\epsilon }\cdot \nabla \phi +{\int }_{\Omega }(1-{\chi }_{\epsilon })(f(u_{\epsilon },w_{\epsilon })-f(u,w))\phi +{\int }_{\Omega }{\partial }_t{W}_{\epsilon }\psi$$

$$+{\int }_{\Omega }(g(u_{\epsilon },w_{\epsilon })-g(u,w))\psi ={\int }_{{\omega }_{\epsilon }}(K_0-K_1)\nabla u\cdot \nabla \phi +{\int }_{{\omega }_{\epsilon }}f(u,w)\phi .$$

$$\frac{1}{2}\frac{d}{dt}\left({\Vert U_{\epsilon }\Vert }_{L^2(\Omega )}^2+{\Vert W_{\epsilon }\Vert }_{L^2(\Omega )}^2\right)+k_{min}{\Vert \nabla U_{\epsilon }\Vert }_{L^2(\Omega )}^2\le 2L\left({\Vert U_{\epsilon }\Vert }_{L^2(\Omega )}^2+{\Vert W_{\epsilon }\Vert }_{L^2(\Omega )}^2\right)$$

$$+{\int }_{{\omega }_{\epsilon }}(K_0-K_1)\nabla u\cdot \nabla U_{\epsilon }+{\int }_{{\omega }_{\epsilon }}f(u,w)U_{\epsilon }.$$

$$\frac{1}{2}\frac{d}{dt}\left({\Vert U_{\epsilon }\Vert }_{L^2(\Omega )}^2+{\Vert W_{\epsilon }\Vert }_{L^2(\Omega )}^2\right)$$

$$\le C\left({\Vert U_{\epsilon }\Vert }_{L^2(\Omega )}^2+{\Vert W_{\epsilon }\Vert }_{L^2(\Omega )}^2\right)+C|{\omega }_{\epsilon }|.$$

Full text

1 Introduction

Cardiac ischemia consists in a restriction of blood supply to the heart tissue usually caused by atherosclerosis or coronary syndrome. The shortage of oxygen may lead to dysfunction of the cell metabolism and eventually to their death. The possible outcomes range from ventricular arrhythmia, fibrillation and ultimately to myocardial infarction.

The ischemic heart syndrome is the most common cardiovascular disease, and the most common cause of death. Hence, the detection of ischemic regions at early stage of their development is of primary importance. This is usually performed by imaging techniques such as echocardiography, gamma ray scintigraphy or magnetic resonance imaging. Nevertheless, the most common test for patients not exhibiting evident symptoms is the electrocardiogram (ECG), which consists in recording electrical impulses across the thorax by means of a set of electrodes. Physicians are often able to identify myocardial ischemia by analysing the evolution of the voltage recorded in the ECG leads although with several technical difficulties (see the mathematical studies in [29] and [10]).

The approach we propose in this paper is to obtain information regarding the electrical functioning of the tissue (and ultimately the presence of ischemic regions) from the knowledge of the electrical potential on the surface of the heart. When employing the potential on the epicardium, i.e., the external boundary of the heart, this can be considered as a natural extension of the well-known mathematical problem often referred to as the inverse problem of electrocardiography, which consists in using the ECG recordings to compute the potential distribution on the epicardium. Alternatively, we remark that it is also possible to obtain electrical measurements on the endocardium, namely the internal boundaries of the heart cavities. Although much more invasive than the ECG, intracardiac ECG (iECG) has become a standard of care in patients with symptoms of heart failure, and allows to get a map of the endocardiac potential by means of non-contact electrodes carried by a catheter inside a heart cavity. We therefore assume that the distribution of the electrical potential is available on some portion of the heart boundary (on the epicardium, in the case of ECG data, or on the endocardium, in the case of iECG measurements). After a reliable model is introduced for the description of the evolution of the electrical potential within the heart, the problem of detecting ischemic regions is hence formulated as an inverse boundary value problem of identifying parameters in a nonlinear reaction diffusion system from boundary data.

In this paper we focus on the determination of small ischemias from boundary measurements of the electrical potential, generalizing the results obtained in previous papers (see [8],[6],[9]) to the monodomain time-dependent model of cardiac electrophysiology.

In general, this model is described by a system consisting in a semilinear parabolic equation coupled with a system of nonlinear ordinary differential equations, where the state variables are the transmembrane potential across the cell membrane, the concentrations of ionic species and the gating variables describing the activity of ionic channels in the membrane. It is well known (see [15] and [11]) that the monodomain model can be derived from the more complex bidomain one, introduced for the first time in [36], for instance by assuming proportionality between intracellular and extracellular conductivity tensors. In particular, in this case, the corresponding conductivity tensor becomes the harmonic mean of the extra and intracellular conductivities. On the other hand, this choice of the conductivity tensor turns out to be the best one to approximate the electrical propagation described by the bidomain model (see, e.g., [28],[16]).

As a first attempt, we limit ourselves to analyze a coupled system of two equations in the potential $u$ and the recovery variable $w$. This corresponds to a large class of phenomenological models, which are characterized by the choice of the nonlinear terms appearing in the partial differential equation and in the ordinary differential equation. Among the most important two-equation systems (see [15] for a general overview) we mention FitzHugh-Nagumo, Rogers-McCulloch and Aliev-Panfilov models. Throughout this paper, we focus on the commonly used version of the Aliev-Panfilov model (originally introduced in [1]), even though the analysis could be extended also to the other two.

In order to detect ischemic regions, we extend the approach of [9], determining a rigorous asymptotic expansion for the perturbed boundary potential due to small conductivity anomalies. To accomplish this task, we need an accurate analysis of the well-posedness of the direct problem for the coupled system in the case of discontinuous anisotropic coefficients and suitable regularity estimates for solutions. In particular, we establish a comparison principle for this class of systems, which to our knowledge was not present in the literature. Here, we consider the case of an insulated heart and we assume to have measurements of the potential on a portion of the boundary.

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 Electrical Impedance Tomography (see, e.g., [3],[14],[5]). A similar approach has also been used in Thermal Imaging (see, e.g., [4]). Here, we are able to extend in a non trivial way the results obtained previously in [6] and [9] for simplified versions of the monodomain model making use of fine regularity for the solutions to nonlinear reaction diffusion systems. In particular, we establish a rigorous expansion for the perturbed transmembrane potential and use this to implement an effective reconstruction algorithm.

The paper is organized as follows. In Section 2 we analyze the well-posedness of the forward problem in the unperturbed and perturbed case. Section 3 is devoted to obtaining energy estimates on the difference between the perturbed and unperturbed electrical potential and an asymptotic expansion of suitable integral terms involving such difference on the boundary. In Section 4 we describe our topological-based optimization algorithm and derive rigorously the topological gradient of a suitable mismatch functional. This is obtained by using the results of Section 3 and some interior regularity results for the solution of parabolic systems. Finally, in Section 5 we outline the numerical implementation of the proposed algorithm, relying on the Finite Element Method for the discrete formulation of the problem. A significant set of numerical experiments is provided in order to assess the effectiveness of the reconstruction, even in presence of data noise.

2 Analysis of the direct problem

The well-posedness analysis of the monodomain system has been object of several studies. We refer to [15, Chapter 3] for a general overview. In [12] a result of existence and uniqueness of weak solution is proved for the FitzHugh-Nagumo, the Aliev-Panfilov and the Rogers-MacCulloch models by means of a Faedo-Galerkin procedure. A result of existence of strong solutions, local in time, is also derived. In [37], instead, results of well-posedness are obtained for a wider range of models, on the base of a fixed point argument.

Regarding the regularity of the solutions of the monodomain system, we report a result in [17] for FitzHugh-Nagumo, Aliev-Panfilov and Rogers-MacCulloch models: if no discontinuities are present in the coefficients of the system, existence and uniqueness of strong solutions is guaranteed, locally in time (see for instance [34] and [20]). A comparison principle is also provided, by means of the tool of invariant sets, allowing to prove existence of global solutions. We also report a result of local existence of classical solutions for the bidomain model, recently obtained in [19].

In this section we focus on the monodomain system in the case of smooth diffusion coefficient and reaction term, corresponding to the case of the healthy tissue (unperturbed case) and in the case of discontinuities in the diffusion coefficient and in the reaction term, corresponding to the presence of an ischemia in the heart tissue (perturbed case). We state an existence, uniqueness and comparison result for classical solutions of the unperturbed case (a proof of which, alternative to the approach of [17], is proposed in [31, Chapter 6]) and prove a result regarding existence, uniqueness and regularity of weak solutions in the perturbed case.

In particular, the initial and boundary value problem associated to the unperturbed monodomain system is the following one

[TABLE]

where $\Omega$ is the region occupied by the hearth tissue, $\nu$ is the outward unit normal vector to the boundary $\partial\Omega$, $Q_{T}:=\Omega\times(0,T)$ and $\Gamma_{T}:\partial\Omega\times(0,T)$. A slightly different formulation of the monodomain model (see, e.g., [15]) involves the presence of a source term in the right-hand side of the first equation in (2.1), representing an applied current, during a limited time window (the initial activation of the tissue). We replace the effect of such a current with the presence of a non-null initial value $u_{0}$. Since the modeling differences between the two formulations are negligible after the first instants, we choose the one proposed in (2.1) more suited for the mathematical analysis of the problem.

In presence of an ischemia $\omega\subset\Omega$, the perturbed case is described by the model

[TABLE]

where $\chi_{\omega}$ is the indicator function of $\omega$ and $K_{\omega}=K_{0}-(K_{0}-K_{1})\chi_{\omega}$. We now specify the requirements on the domain, the coefficients and the source terms.

Assumption 1**.**

$\Omega\subset\mathbb{R}^{d}$* bounded domain, $d=2,3$, and $\partial\Omega\in C^{2+\alpha}$;* 2. 2.

the inclusion $\omega\subset\Omega$ is well separated from the boundary, i.e.,

[TABLE] 3. 3.

$K_{0},K_{1}\in C^{2}(\overline{\Omega};\mathbb{R}^{d\times d})$* are symmetric matrix-valued functions in $\Omega$; $\forall x\in\overline{\Omega}$, the matrices $K_{0}(x)$ and $K_{1}(x)$ admit $d$ positive eigenvalues $k_{0,1}\leq\ldots\leq k_{0,d}$ and $k_{1,1}\leq\ldots\leq k_{1,d}$ respectively, associated to the same eigenvectors $\vec{e}_{1}(x),\ldots\vec{e}_{d}(x)$ such that $k_{1,i}\leq k_{0,i}$ $\forall i=1,\ldots,d$;* 4. 4.

$u_{0}\in C^{2+\alpha}(\overline{\Omega})$, $w_{0}\in C^{2}(\overline{\Omega})$, and $K_{0}\nabla u_{0}\cdot\nu=0$ on $\partial\Omega$; 5. 5.

we assume $f,g$ as in the Aliev-Panfilov model, namely:

[TABLE]

being $A,\epsilon>0$ and $0<a<1$.

Remark 1**.**

The requirements on the matrix-valued functions $K_{0}$ and $K_{1}$ are satisfied by the conductivity tensors prescribed in the model under consideration (see [15, 35]). According to experimental evidence, the cardiac tissue can be modeled as an orthotropic material, characterized by the presence of fibers and sheets, which define the conductivity eigenvectors. Moreover, the presence of an ischemia does not affect the direction of the fibers, but reduces the value of all the associated eigenvalues. From now on, we will indicate by $k_{min}$ the minimum eigenvalue of $K_{1}$ and by $k_{max}$ the maximum eigenvalue of $K_{0}$.

Remark 2**.**

As an immediate consequence of (2.4), $f$ and $g$ satisfy the Tangency condition on the rectangle $S\mathrel{\mathop{:}}=[0,1]\times[0,\frac{A(a+1)^{2}}{4}]$, (see [2]), i.e., indicating by $\vec{p}$ a generalized outward normal on $\partial S$ ( for all $(\xi_{1},\xi_{2})\in\partial S$, $\vec{p}(\xi_{1},\xi_{2})$ is such that $\vec{p}(\xi_{1},\xi_{2})\cdot(\xi_{1},\xi_{2})\geq\vec{p}(\xi_{1},\xi_{2})\cdot(\eta_{1},\eta_{2})$ $\forall(\eta_{1},\eta_{2})\in S$) then,

[TABLE]

Moreover, the functions $f,g$ are Lipschitz continuous on $S$ with constants $L_{f},L_{g}\leq L$.

Remark 3**.**

For the sake of brevity, in all the formulas we avoid to indicate time and space integration variables with respect to the classical Lebesgue measure, unless it is necessary.

We now outline the main results regarding the well-posedness of the problems (2.1) and (2.2).

Theorem 1** (Unperturbed problem).**

Let Assumption 1 holds, and suppose that, $\forall x\in\overline{\Omega},$ $(u_{0}(x),w_{0}(x))\in S$. Then, problem (2.1) admits a unique classical solution $(u,w)$, namely $u\in C^{2+\alpha,1+\alpha/2}(\overline{Q_{T}})$, $w\in C^{\alpha,1+\alpha/2}(\overline{Q_{T}})$. Moreover, $(u(x,t),w(x,t))\in S$, for each $(x,t)\in\overline{Q_{T}}$.

The proof of Theorem 1 is derived by means of classical fixed point argument (see [27, Chapter 8, Sections 9 and 11]). A detailed proof is reported in [31, Theorem 6.1].

Regarding the perturbed problem (2.2), we note that, although the conductivity tensor and the nonlinear term are discontinuous, we can extend the results obtained in [12] thanks to the uniform ellipticity to the boundedness of the conductivity tensors and to the form of the reaction term, deriving the following existence and uniqueness result.

Theorem 2** (Perturbed problem).**

Let the Assumption 1 holds, and suppose that, $\forall x\in\overline{\Omega},$ $(u_{0}(x),w_{0}(x))\in S$. Then, problem (2.2) admits a unique weak solution, i.e., a couple $(u_{\omega},w_{\omega})$ such that $u_{\omega}\in L^{2}(0,T;H^{1}(\Omega))\cap L^{\infty}(0,T;L^{2}(\Omega))$, $\partial_{t}u_{\omega}\in L^{2}(0,T;H^{*})$, $w_{\omega}\in L^{\infty}(0,T;L^{2}(\Omega))$, $\partial_{t}w_{\omega}\in L^{2}(Q_{T})$ and satisfying, for a.e. $t\in(0,T)$,

[TABLE]

Moreover, $u_{\omega}\in C^{\alpha,\alpha/2}(\overline{Q_{T}})$, $w_{\omega}\in C^{\alpha,\alpha/2}(\overline{Q_{T}})$ and $(u_{\omega}(x,t),w_{\omega}(x,t))\in S$, $\forall\,(x,t)\in\overline{Q_{T}}$.

Proof.

To start with, note that uniqueness of the weak solution of (2.2) has been shown in the case of Aliev-Panfilov model in [21, Theorem 1.3] as a byproduct of a stability result obtained by exploiting the specific nonlinear expression of $f$ and $g$.

We proceed now introducing a sequence of regularized problems $(P_{k}),\,k\in\mathbb{N},$ of (2.2) and showing that the sequence of their solutions converges to a weak solution of (2.2). We then exploit additional properties inherited by the approximation process to conclude the stated regularity results. The uniqueness argument is briefly sketched.

Since $\chi_{\omega}$ is an indicator function, surely $\chi_{\omega}\in L^{2}(\Omega)$; by density arguments and according to (2.3),

[TABLE]

being $C^{2}_{c}(\Omega)$ the space of $C^{2}$ functions with compact support in $\Omega$. Define $K_{k}=K_{0}+(K_{1}-K_{0})\phi_{k}$ and let $(u_{k},w_{k})$ be the solution of the following problem

[TABLE]

where we have used the fact that $K_{k}=K_{0}$ on $\partial\Omega$. We observe that, for any fixed $k$, an application of Theorem 1 ensures the existence and uniqueness of a classical solution of Problem (2.8). Moreover by (2.7) we have that $(1-\phi_{k})f$ and $g$ satisfy the Tangency condition on $S$. Hence, from Theorem 1 we deduce $(u_{k},w_{k})\in S$, for all $k$.

We now prove that from $\phi_{k}\xrightarrow{L^{2}}\chi_{\omega}$ the convergence of $(u_{k},w_{k})$ to a weak solution $(u,w)$ of (2.2) holds. We start by proving some a priori estimates. Consider the weak form of the problem solved by $(u_{k},w_{k})$ and take the classical solutions $u_{k}$, $w_{k}$ as test functions. Then, we get

[TABLE]

Recall now that $k_{min}$ is the minimum eigenvalue of $K_{1}$, whereas $k_{max}$ is the maximum eigenvalue of $K_{0}$. Moreover, since $\phi_{k}$, $u_{k}$, $w_{k}$ are uniformly bounded independently of $k$ (indeed, $\phi_{k}(x)\in[0,1]$ and $(u_{k},w_{k})\in S$) and $f,g$ are continuous, we can introduce $M_{f}:=max_{(x,t)\in Q_{T}}\left|(1-\phi_{k})f(u_{k},w_{k})\right|$ and $M_{g}:=max_{(x,t)\in Q_{T}}\left|g(u_{k},w_{k})\right|$, which are independent of $k$. Hence, by Young inequality,

[TABLE]

Using Gronwall’s inequality, we get

[TABLE]

which implies

[TABLE]

Integrating from [math] to $t$ (2.9) and using last estimate it also follows that

[TABLE]

A bound for the $H^{*}$ norm of $\partial_{t}u$ can be found by considering that, for each $\varphi\in H^{1}(\Omega)$,

[TABLE]

and computing the $L^{2}$ norm in time

[TABLE]

Analogously, one proves that ${\left\|\partial_{t}w_{k}\right\|}_{L^{2}(Q_{T})}\leq c_{4}$, with $c_{4}$ independent of $k$.

As a consequence of the uniform bounds we can ensure that (up to a subsequence) $\exists u\in L^{2}(0,T;H^{1}(\Omega))\cap L^{\infty}(0,T;L^{2}(\Omega))$, $\exists w\in L^{\infty}(0,T;L^{2}(\Omega))$, $\exists u^{*}\in L^{2}(0,T;H^{*})$, $\exists w^{*}\in L^{2}(Q_{T})$ such that

[TABLE]

Using the result contained in [33, Theorem 8.1] and the uniqueness of the weak solution of the perturbed problem we have that $u_{k}\xrightarrow{L^{2}(Q_{T})}u_{\omega}$ (see [33, Theorem 8.1]), hence $u_{k}\rightarrow u_{\omega}$ a.e. in $Q_{T}$. Moreover, also $w_{k}\rightarrow w_{\omega}$ a.e. in $Q_{T}$ as it can be straightforwardly obtained by the expression

[TABLE]

Consider now the weak form of the problem solved by $(u_{k},w_{k})$: $\forall\varphi\in H^{1}(\Omega)$, $\psi\in L^{2}(\Omega)$.

[TABLE]

Taking the limit in $k$, exploiting the weak convergence of $u_{k}$, $w_{k}$, $\partial_{t}u_{k}$, $\partial_{t}w_{k}$, we obtain that $u^{*}=\partial_{t}u_{\omega}$, $w^{*}=\partial_{t}w_{\omega}$ in $\mathcal{D}^{\prime}(0,T)$ and that the limit $(u_{\omega},w_{\omega})$ satisfies, $\forall\varphi\in H^{1}(\Omega)$, $\psi\in L^{2}(\Omega)$,

[TABLE]

Indeed, the convergence of the terms involving the time derivatives is a direct consequence of the weak convergence of $\partial_{t}u_{k}$, $\partial_{t}w_{k}$ and of the definition of distributional derivative. The limit of the nonlinear reaction terms can be proved by taking advantage of the dominated convergence theorem and of the pointwise (a.e.) convergence of $u_{k}$ and $\phi_{k}$. The convergence of the diffusion term is obtained by combining the weak convergence of $u_{k}$ in $H^{1}(\Omega)$, the pointwise (a.e.) convergence of $\phi_{k}$ and the uniform $L^{\infty}(\Omega)$ bound on $\phi_{k}$. According to (2.12), we can ensure that the limit $(u_{\omega},w_{\omega})$ is the weak solution of (2.2).

The weak solution $(u_{\omega},w_{\omega})$ is moreover a pointwise (a.e.) limit of the regularized solutions $(u_{k},w_{k})$. As a consequence, the uniform bound on $(u_{k},w_{k})$ is valid also for the limit: $(u_{\omega}(x,t),w_{\omega}(x,t))\in S$ a.e. in $Q_{T}$. This allows to prove the additional Hölder regularity of $u_{\omega}$. Indeed, since $K_{\omega}\in L^{\infty}(\Omega)$, $f(u_{\omega},w_{\omega})\in L^{\infty}(Q_{T})$, we can apply Theorem 10.1 of [22, Chapter 3] on the first equation in (2.2)

[TABLE]

to get $u_{\omega}\in L^{\infty}(Q_{T})$. Now, we can extend the regularity result up to the boundary due to the hypothesis on $\partial\Omega$ and on $u_{0}$ contained in Assumption 1, and conclude $u_{\omega}\in C^{\alpha,\alpha/2}(\overline{Q_{T}})$. Using the analytic expression of $w_{\omega}$ that can be obtained by (2.10) and the regularity of $u_{\omega}$, we can also deduce that $w_{\omega}\in C^{\alpha,1+\alpha/2}(\overline{Q_{T}})$. ∎

3 Analysis of the inverse problem

We now tackle the inverse problem of identifying a perturbation $\omega$ from boundary measurements. Suppose to know $u_{meas}$, the trace of the solution of (2.2) in presence of an unknown inclusion. The inverse problem reads

[TABLE]

Although the analysis of the direct problem has been performed in presence of an arbitrary inclusion $\omega\subset\Omega$, for the purpose of solving the inverse problem and derive a reconstruction algorithm, we limit ourselves at considering the case of inclusions of small size, in analogy to what done in [8], [9]. In particular, we consider a family of inclusions $\omega_{\varepsilon}$ satisfying (2.3) for each $\varepsilon$ and such that

[TABLE]

We define $\chi_{\varepsilon}$ the indicator function of $\omega_{\varepsilon}$, $K_{\varepsilon}:=K_{0}-(K_{0}-K_{1})\chi_{\varepsilon}$ and $(u_{\varepsilon},w_{\varepsilon})$ the solution of problem (2.2) with $\omega=\omega_{\varepsilon}$, in the sense of Theorem 2.

3.1 Energy estimates

We now derive some energy estimates for the difference between $(u_{\varepsilon},w_{\varepsilon})$ and $(u,w)$, the solution of the unperturbed problem (2.1) in terms of $|\omega_{\varepsilon}|$ when $\varepsilon\rightarrow 0$.

From now on we will indicate by $C$ a positive constant depending on the data, independent of $\varepsilon$ and that may vary also in the same line.

Proposition 1**.**

Under Assumption 2.1 the following inequalities hold

$\displaystyle{\left\|u_{\varepsilon}-u\right\|}_{L^{\infty}(0,T,L^{2}(\Omega))}+{\left\|w_{\varepsilon}-w\right\|}_{L^{\infty}(0,T,L^{2}(\Omega))}\leq C|\omega_{\varepsilon}|^{\frac{1}{2}}$, $\displaystyle{\left\|u_{\varepsilon}-u\right\|}_{L^{2}(0,T,H^{1}(\Omega))}\leq C|\omega_{\varepsilon}|^{\frac{1}{2}}$, $\displaystyle{\left\|u_{\varepsilon}-u\right\|}_{L^{2}(Q_{T})}\leq C|\omega_{\varepsilon}|^{\frac{1}{2}+\eta},\quad\text{ for some }\eta>0$.

Proof.

We consider (2.6) with $\omega=\omega_{\varepsilon}$ and the weak formulation of the unperturbed problem (2.1). Subtracting term by term and defining ${U_{\varepsilon}}=u_{\varepsilon}-u$ and ${W_{\varepsilon}}=w_{\varepsilon}-w$, we have that ${U_{\varepsilon}}(\cdot,0)=0$, ${W_{\varepsilon}}(\cdot,0)=0$ and, for almost every $t\in(0,T)$,

[TABLE]

Let $\varphi={U_{\varepsilon}}$, $\psi={W_{\varepsilon}}$ then, according to Theorem 1 and Theorem 2, both $(u,w)$ and $(u_{\varepsilon},w_{\varepsilon})$ range within the rectangle $S$, on which the functions $f$ and $g$ are Lipschitz continuous with constants less or equal than $L$. Hence

[TABLE]

Via an application of Schwarz and Young inequalities on the last two terms in the right-hand side of the previous inequality, and from the regularity of the solution $(u,w)$, we can deduce

[TABLE]

An application of Gronwall’s lemma to (3.4) entails that ${\left\|{U_{\varepsilon}}\right\|}_{L^{\infty}(0,T;L^{2}(\Omega))}^{2}$, ${\left\|{W_{\varepsilon}}\right\|}_{L^{\infty}(0,T;L^{2}(\Omega))}^{2}$ and ultimately ${\left\|{U_{\varepsilon}}\right\|}_{L^{2}(0,T;H^{1}(\Omega))}^{2}$ can be bounded by $C|\omega_{\varepsilon}|$. This allows to conclude the first two statements.

Observe now that the pair $({U_{\varepsilon}},{W_{\varepsilon}})$ is also the solution of

[TABLE]

By the mean value theorem, there exist two pairs of functions $(u_{\xi_{1}},w_{\xi_{1}})$, $(u_{\xi_{2}},w_{\xi_{2}})$ s.t. ***it follows by an application of the Lagrange’s mean value theorem on the real-valued function $h(\tau)=\int_{\Omega}(1-\chi_{\varepsilon})(f(u+\tau(u_{\varepsilon}-u),w+\tau(w_{\varepsilon}-w))-f(u,w)){U_{\varepsilon}}$ on the interval $[0,1]$; the same holds for $g$.

[TABLE]

By definition, $(u_{\xi_{1}},w_{\xi_{1}})$ and $(u_{\xi_{2}},w_{\xi_{2}})$ are convex combinations of $(u,w)$ and $(u_{\varepsilon},w_{\varepsilon})$, thus they assume values in the rectangle $S$. Let now $(\overline{U}_{\varepsilon},\overline{W}_{\varepsilon})$ be the solution of the adjoint problem

[TABLE]

with initial conditions $\overline{U}_{\varepsilon}(\cdot,T)=0$, $\overline{W}_{\varepsilon}(\cdot,T)=0$ and homogeneous Neumann boundary condition. Consider the change of variable: $z(\cdot,t)=\overline{U}_{\varepsilon}(\cdot,T-t)$, $y(\cdot,t)=\overline{W}_{\varepsilon}(\cdot,T-t)$ and define $\widehat{U}_{\varepsilon}(\cdot,t)={U_{\varepsilon}}(\cdot,T-t)$, $\widehat{W}_{\varepsilon}(\cdot,t)={W_{\varepsilon}}(\cdot,T-t)$, $\hat{f}_{u}(\cdot,t)=f_{u}(u_{\xi_{1}}(\cdot,T-t),w_{\xi_{1}}(\cdot,T-t))$, and analogously for $\hat{f}_{w},\hat{g}_{u},\hat{g}_{w}$. Hence $z$ and $y$ solve

[TABLE]

Since $\widehat{U}_{\varepsilon},\widehat{W}_{\varepsilon}\in L^{2}(Q_{T})$ and $\hat{f}_{u},\hat{f}_{w},\hat{g}_{u},\hat{g}_{w}$ are bounded in $\overline{Q}_{T}$, by standard Faedo-Galerkin technique we obtain that the solution $(z,y)$ of (3.7) exists and is unique with the properties $z\in L^{2}(0,T;H^{1}(\Omega))\cap L^{\infty}(0,T;L^{2}(\Omega))$, $\partial_{t}z\in L^{2}(0,T;H^{*})$, $y\in L^{\infty}(0,T;L^{2}(\Omega))$, $\partial_{t}y\in L^{2}(Q_{T})$. Moreover

[TABLE]

Additional regularity of $z$ can be proved with an analogous argument as in [22, Chapter 4, Theorem 9.1], applied on the first equation in (3.7). Indeed, by the regularity of $K_{0}$ and the square integrability of $\widehat{U}_{\varepsilon}-\hat{g}_{u}y$, we can conclude that $z\in L^{2}(0,T;H^{2}(\Omega))$, $\partial_{t}z\in L^{2}(Q_{T})$ and

[TABLE]

Moreover, multiplying the first two equations in (3.7) respectively by $\partial_{t}z$ and $\partial_{t}y$ and integrating on $\Omega$, straightforward computations allow to conclude that $z\in L^{\infty}(0,T;H^{1}(\Omega))$ with

[TABLE]

By Sobolev inequality we obtain

[TABLE]

(all the norms being bounded by ${\left\|{U_{\varepsilon}}\right\|}_{L^{2}(Q_{T})}+{\left\|{W_{\varepsilon}}\right\|}_{L^{2}(Q_{T})}$; the same bounds hold on $\overline{U}_{\varepsilon}$). Thus, via Nirenberg interpolation estimates (see [25]), we get

[TABLE]

Recalling also the previous results, this allows to conclude that, taking $p\in\left(2,\frac{10}{3}\right]$,

[TABLE]

Let us multiply the equations of (3.5) respectively by $\overline{U}_{\varepsilon},\overline{W}_{\varepsilon}$ and the first two equations of (3.7) respectively by ${U_{\varepsilon}},{W_{\varepsilon}}$. Integrating on $(0,T)$ and summing the resulting identities it is easy to see that

[TABLE]

Thanks to (3.8) and Hölder inequality the first term of the right-hand side can be bounded as follows

[TABLE]

where $p\in\left(2,\frac{10}{3}\right]$ and $q=\frac{p}{p-1}\in\left[\frac{10}{7},2\right)$. In addition, again by Hölder inequality,

[TABLE]

Furthermore, it is straightforward to see that the last term in (3.9) can be bounded by $c|\omega_{\varepsilon}|^{\frac{1}{q}}{\left\|\overline{U}_{\varepsilon}\right\|}_{L^{2}(Q_{T})}$. Finally, from (3.9) we conclude that, for $q\in\left[\frac{10}{7},2\right)$,

[TABLE]

Hence

[TABLE]

where $\beta\in\left(0,\frac{1}{5}\right]$. ∎

3.2 Asymptotic expansion of boundary voltage

In this section we prove the following result.

Theorem 3**.**

For every sequence $\{\omega_{\varepsilon_{n}}\}$ with $|\omega_{\varepsilon_{n}}|>0$, $|\omega_{\varepsilon_{n}}|\rightarrow 0$, there exist a subsequence (still denoted by $\omega_{\varepsilon_{n}}$), a Radon measure $\mu$ and a matrix-valued function $\mathcal{M}\in L^{2}(\Omega,d\mu;\mathbb{R}^{d\times d})$ such that, for every pair $(\Phi,\Psi)\in(C^{1}(\overline{Q}_{T}),C(\overline{Q}_{T}))$ satisfying

[TABLE]

with the final conditions $\Phi(\cdot,T)=0$, $\Psi(\cdot,T)=0$, it holds

[TABLE]

In order to prove Theorem 3.11, we need to introduce the auxiliary functions ${v^{(j)}}$, ${v^{(j)}_{\varepsilon_{n}}}$ solving

[TABLE]

with $F^{(j)}=\sum_{i}\frac{\partial[K_{0}]_{ij}}{\partial x_{i}}$ and $f^{(j)}=\sum_{i}[K_{0}]_{ij}\nu_{i}$. By the choice of $F^{(j)}$ and $f^{(j)}$ it holds ${v^{(j)}}=x_{j}-\int_{\partial\Omega}x_{j}$. Moreover, ${v^{(j)}}$ and ${v^{(j)}_{\varepsilon_{n}}}$ satisfy the following energy estimates (see [13])

[TABLE]

We need also a preliminary lemma.

Lemma 1**.**

For each $\phi\in C^{1}(\overline{Q}_{T})$ such that $\phi(x,T)=0$, we have, as $|{\omega_{\varepsilon_{n}}}|\rightarrow 0$,

[TABLE]

Proof.

Since ${v^{(j)}}$ and ${v^{(j)}_{\varepsilon_{n}}}$ are the solutions of problems (3.12), we obtain

[TABLE]

Take $\varphi={U_{\varepsilon_{n}}}\phi$, being ${U_{\varepsilon_{n}}}={u_{\varepsilon_{n}}}-u$. By computation we get

[TABLE]

that we can rewrite as

[TABLE]

Proposition 1 and (3.13) lead to

[TABLE]

whereas,

[TABLE]

Collecting the previous two relations in (3.15), we conclude

[TABLE]

It can be easily verified that, for every $\varphi\in H^{1}(\Omega)$, the following identities hold

[TABLE]

Taking $\varphi={v^{(j)}_{\varepsilon_{n}}}\phi$ in the first identity and $\varphi={v^{(j)}}\phi$ in the second one, integrating in time and subtracting we obtain

[TABLE]

Via integration by parts, and taking advantage of the homogeneous initial and final conditions satisfied respectively by $U_{\varepsilon_{n}}$ and $\phi$ it follows

[TABLE]

Analogous bounds can be proved for the last term in the left-hand side and the first term in the right-hand side of (3.17), thanks to the energy estimates of ${v^{(j)}_{\varepsilon_{n}}}-{v^{(j)}},\,{u_{\varepsilon_{n}}}-u,\,{w_{\varepsilon_{n}}}-w$ and the regularity of $f$. As a consequence, it holds that

[TABLE]

In conclusion, by a combination of (3.16) and (3.18),

[TABLE]

which immediately entails the thesis. ∎

We are finally ready to prove the main result of this section.

Proof of Theorem 3.11.

Arguing as in [13], we deduce that there exist a Radon measure $\mu\in C(\overline{\Omega})^{*}$, a symmetric matrix $\mathcal{M}\in L^{2}(\Omega,d\mu;\mathbb{R}^{d\times d})$ and a sequence $\{{\omega_{\varepsilon_{n}}}\}$ such that

[TABLE]

in the weak* topology of $C(\overline{\Omega})$. Moreover, thanks to the regularity estimates for $u$ and the symmetries of the matrices $K_{0}$, $K_{1}$ and $\mathcal{M}$, we have, for every $j=1,\ldots,d$,

[TABLE]

which also implies

[TABLE]

in the weak∗ topology of $C(\overline{Q}_{T})$.

On account of the energy estimates (3.13), straightforward computations show that for every $i,j=1,\ldots,d$ there exists a constant $C$, depending on the data and on $T$ but independent of $\varepsilon_{n}$, such that

[TABLE]

In particular, we define the limit measure $\nu_{j}$ in the weak∗ topology of $C(\overline{Q}_{T})$) as follows

[TABLE]

Exploiting (3.21) and (3.20), and recalling Lemma 3.3, we deduce

[TABLE]

and by the density of $C^{1}(\overline{Q}_{T})$ in $C^{0}(\overline{Q}_{T})$ we derive

[TABLE]

Now, setting ${U_{\varepsilon_{n}}}={u_{\varepsilon_{n}}}-u$ and ${W_{\varepsilon_{n}}}={w_{\varepsilon_{n}}}-w$, and selecting $\Phi$ and $\Psi$ as test functions, we have

[TABLE]

Furthermore, since $(\Phi,\Psi)$ is a solution to (3.10), it also holds:

[TABLE]

Summing up the last two identities and integrating in time we observe that the terms involving the time derivatives vanish due to the initial and final conditions on ${U_{\varepsilon_{n}}},{W_{\varepsilon_{n}}},\Phi,\Psi$. Consider now the non-linear terms containing $f$

[TABLE]

Observe that by means of the Lagrange’s mean value theorem and the Lipschitz-continuity of the functions $f_{u},f_{w}$ with constants $L_{f_{u}},L_{f_{w}}$ we have

[TABLE]

whereas

[TABLE]

Analogous estimates can be obtained for the terms containing $g$. Finally,

[TABLE]

Thanks to the regularity of $\Phi$, and employing the computed weak* limits, we derive

[TABLE]

which concludes the proof. ∎

4 Reconstruction algorithm: a topology optimization approach

The asymptotic expansion provided in Theorem 3.11 allows to describe the perturbation of the electrical potential on the boundary of the domain due to the presence of a small conductivity inhomogeneity $\omega_{\varepsilon}$. In order to derive a reconstruction algorithm for problem (3.1), we introduce the mismatch functional $J$

[TABLE]

where $u_{\varepsilon}=u_{\omega_{\varepsilon}}$ solves the perturbed problem (2.2) in the presence of an ischemic region $\omega_{\varepsilon}$. We now prove that the functional $J$ restricted to the class of inclusions satisfying (2.3), (3.2) admits an asymptotic expansion with respect to the size of the inclusion. Moreover, as it is shown in Theorem 4.3, the first-order term of the expansion (which will be denoted as $G$, the topological gradient of $J$) can be computed by solving the unperturbed problem and a suitable adjoint problem. In particular, we restrict ourselves to inclusions $\omega_{\varepsilon}$ satisfying (2.3) and of the form

[TABLE]

where $B$ is a bounded, smooth set containing the origin.

We have the following

Theorem 4**.**

Consider a family $\{\omega_{\varepsilon}\}$ satisfying (2.3) and (4.2). Then, there exists a matrix $\mathcal{M}$ (which may depend on $z$, $B$, $K_{0}$ and $K_{1}$) such that, as $\varepsilon\to 0$,

[TABLE]

where $(u,w)$ solves (2.1) and $(\Phi,\Psi)$ is the solution of the adjoint problem:

[TABLE]

Proof.

First of all, we note that, since Lemma 3.3 holds under the weaker assumption on the test functions, $\phi\in W^{1,\infty}({Q}_{T})$, then the validity of Theorem 3.11 can be easily extended when $(\Phi,\Psi)\in(W^{1,\infty}({Q}_{T}),C(\overline{Q}_{T}))$. We now show that the solution to (4.3) enjoys this regularity. In fact this can be ensured by a lifting argument, since the boundary datum $u-u_{meas}\in C^{2+\alpha,1+\alpha/2}(\Gamma_{T})$. Estimates of the $W^{1,\infty}(Q_{T})$ norm of $\Phi$ can be obtained by flattening the boundary and by reflection arguments, see for example [22, Chapter 7, Theorem 2.2] for a two-dimensional case.

It is straightforward to verify that

[TABLE]

The first term in the right-hand side of (4.4) involves the boundary condition of the adjoint problem, and in particular it can be written as $\int_{0}^{T}\int_{\partial\Omega}K_{0}\nabla\Phi\cdot\nu(u_{\varepsilon}-u)$. We now apply Theorem 3.11 in order to conclude that

[TABLE]

According to assumption (4.2), the limit measure $\mu$ is independent of the choice of the subsequence and is equal to $\delta_{z}$, the Dirac measure centered in $z$. Finally, the second term in the right-hand side of (4.4) is $o(|\omega_{\varepsilon}|)$ by means of Lemma 2 below. ∎

Remark 4**.**

By Assumption 2.1 and by the homogeneous boundary conditions on $u_{\varepsilon}$ we can conclude, using standard local regularity results on parabolic equations applied in a neighbourhood of $\partial\Omega$, that $u_{\varepsilon}\in C^{2+\alpha,1+\alpha/2}(\Gamma_{T})$

Lemma 2**.**

Let $\{\omega_{\varepsilon}\}$ satisfy (2.3) and (4.2). Then,

[TABLE]

Proof.

By the same argument as in the proof of Theorem 4.3 and by Remark 4.2,

[TABLE]

where $(\Theta,\Xi)$ satisfies

[TABLE]

and $(\Theta,\Xi)\in(W^{1,\infty}({Q}_{T}),C(\overline{Q}_{T}))$.

We now focus on proving that for a compact set $K$ such that $\omega_{\varepsilon}\subset K$ and satisfying $dist(K,\partial\Omega)\geq d_{0}$ one has

[TABLE]

and, as a consequence of a Sobolev immersion, (4.6) implies (here $d=2,3$)

[TABLE]

By a combination of (2.3) and (4.2) we can ensure that the support of the measure $\mu=\delta_{z}$ is contained in $K$. This would entail that

[TABLE]

hence $\int_{0}^{T}\int_{\partial\Omega}(u_{\varepsilon}-u)^{2}=O(|\omega_{\varepsilon}|^{3/2})=o(|\omega_{\varepsilon}|)$.

In order to prove (4.6), by standard Faedo-Galerkin technique applied to the linear system (4.5), one has

[TABLE]

Moreover, according to (4.5), $\Theta$ is the solution of the following problem

[TABLE]

Analogously to what stated in [22, Chapter 4, Theorem 9.1], it holds that

[TABLE]

and in conclusion, by the regularity of $g,u,w$ and by (4.7)

[TABLE]

Similar estimates can be derived also for $\Xi$ since $f,g\in C^{3}(\mathbb{R}^{2})$ and $u,w\in C^{2,1}(\overline{Q_{T}})$. This last property is guaranteed on $u$ by Theorem 1, whereas it can be extended to $w$ since we consider $w_{0}\in C^{2}(\overline{\Omega})$. Now, by computing the solution of the third equation in (4.5) in closed form, we can easily verify that

[TABLE]

Consider the change of variable $s=T-t$ and denote by $\widehat{\Theta}(\cdot,s)=\Theta(\cdot,T-t),\widehat{\Xi}(\cdot,s)=\Xi(\cdot,T-t)$. We now focus on the first equation in (4.5) and compute the second derivatives of each term. Then, $V\mathrel{\mathop{:}}=\frac{\partial^{2}\widehat{\Theta}}{\partial x_{i}\partial x_{j}}$ satisfies the following equation in a weak sense:

[TABLE]

being

[TABLE]

Let $\Omega_{1}$ be an open subset of $\Omega$, such that $\omega_{\varepsilon}\subset K\subset\Omega_{1}\subset\subset\Omega$. By interior regularity results it holds that $(\widehat{\Theta},\widehat{\Xi})$ is smooth on $\Omega_{1}\times[0,T]$, which entails by the initial conditions that

[TABLE]

According to the smoothness of $f,g,u,w$ and to the $H^{2}$ bounds (4.8) and (4.9), we get

[TABLE]

In fact, all the terms in $R$ except the first three belong to $L^{2}(Q_{T})$, with norm bounded by ${\left\|\Theta\right\|}_{L^{2}(0,T;H^{2}(\Omega))}$ or by ${\left\|\Xi\right\|}_{L^{2}(0,T;H^{2}(\Omega))}$. Set $H:=div\left(\frac{\partial K_{0}}{\partial x_{i}}\nabla\frac{\partial\widehat{\Theta}}{\partial x_{j}}\right)$, then, for any $v\in H^{1}_{0}(\Omega_{1})$,

[TABLE]

where $\langle\cdot,\cdot\rangle_{\star}$ indicates the duality pairing between the involved spaces. Integrating in time last relation we finally get

[TABLE]

With similar arguments we can estimate also the second and third term in the definition of the function $R$ leading finally to (4.12).

Consider now a test function $\zeta\in C^{\infty}_{C}(\Omega_{1})$, $0\leq\zeta\leq 1$, such that $\zeta=1$ in $K$. According to (4.10), by simple computations we verify that $\widetilde{V}=V\zeta$ satisfies

[TABLE]

being $\widetilde{R}=\zeta R-2div(VK_{0}\nabla\zeta)+Vdiv(K_{0}\nabla\zeta)$. It holds that $div(VK_{0}\nabla\zeta)\in L^{2}(0,T;H^{-1}(\Omega_{1}))$, and

[TABLE]

hence we observe that

[TABLE]

By standard Faedo-Galerkin argument (see, e.g., [18, Chap. 18, Par. 3, Theorem 3], we verify that

[TABLE]

Finally, by the definition of $\zeta$, combining (4.15) and (4.14) we conclude that

[TABLE]

The thesis is now proved thanks to the trace inequality and the energy estimates. ∎

Theorem 4.3 gives a representation formula for the topological gradient $G$, the first order term in the expansion of the mismatch functional $J$. In analogy to [6, 9], we propose a one-step reconstruction Algorithm 1 for the identification of small inclusions satisfying (4.2).

Remark 5**.**

The polarization tensor $\mathcal{M}$ can be computed in an explicit form, e.g., when the shape $B$ of the inclusion is a disk, even though it is not straightforward in the current anisotropic case. Since $\mathcal{M}$ depends on $K_{0}$ and $K_{1}$ and on the shape of the set $B$, which we assume to be a disk, it varies according to the space variable $z$. In order to compute $\mathcal{M}(z)=\mathcal{M}(K_{0}(z),K_{1}(z),B)$, we first use [3, Lemma 4.30] to transform $K_{0}$ and $K_{1}$ in diagonal matrices. By a rescaling of the spatial variables, we can reduce to the case in which one of the two diagonalized tensors is the identity, which allows to apply [3, Proposition 4.31]. We remark that the anisotropic rescaling entails that the original circle $B$ is transformed in an ellipse of known semiaxes.

5 Numerical results

In this section we describe the implementation of Algorithm 1 and report some numerical results in order to show its effectiveness. In particular, we set our experiments in a two-dimensional idealized geometry, representing a horizontal section of the ventricles. The application of the model on a three-dimensional geometry would be equivalently possible, and will be object of future studies, employing advanced numerical analysis techniques in order to tackle the increased computational effort. We rely on synthetic data, i.e., we solve the monodomain system in presence of a prescribed ischemic region and then use the value of the solution on the boundary (or a portion of it) as an input for the reconstruction algorithm. In order to prevent inverse crimes, we employ different numerical settings for the synthetic data generation and for the solution of the unperturbed and adjoint problems, required for the reconstruction algorithm. In particular, we adopt a much more refined discretization, both in space and in time, for the simulation of the synthetic data.

In the following experiments we assume to measure the voltage only on a portion of the heart. As outlined in Section 1, measurements of the voltage can be acquired on the inner surface of a ventricle by intracavitary measurements, or, alternatively, we might be able to compute a map of the electrical potential on the epicardium starting from ECG data. This does not affect the reconstruction procedure described in Algorithm 1, apart from the definition of the adjoint problem (4.3), which now prescribes oblique boundary conditions involving $u_{meas}-u$ on the portion of the boundary where the measurements are acquired, and homogeneous Neumann conditions elsewhere.

5.1 Finite Element approximation

In order to numerically approximate the solution of the monodomain model, we rely on a Galerkin Finite Element scheme, introducing a tessellation of the domain consisting of triangular elements. In particular, we adopt two different meshes for the solution of the unperturbed (and adjoint) problem and for the generation of the synthetic measurements, the latter being more refined especially close to the boundary of the prescribed ischemic region.

Moreover, in order to reproduce the anisotropic behavior of the conductivity coefficients $K_{0}$ and $K_{1}$, we consider the presence of fibers within the domain. We adopt a procedure analogous to the one reported in [30] for the generation of the fiber directions, resorting on the solution of a Laplace problem with suitable boundary conditions. Once the direction of the fibers $\vec{e_{f}}(x)=\vec{e}_{1}(x)$ is defined within $\Omega$, as well as the transmural vector field $\vec{e_{n}}(x)=\vec{e}_{2}(x)$ (obtained by a clockwise rotation of $90^{\circ}$ of $\vec{e_{f}}(x)$), the definition of the conductivity tensors is

[TABLE]

The numerical mesh for the solution of the background problem and the directions of the fibers are reported in Figure 1.

As already mentioned in the introduction, the best choice of the conductivity tensor $K_{0}$ to be employed in the numerical simulations of the monodomain model is given by the harmonic mean of the intracellular and extracellular conductivity tensors $D_{i}$ and $D_{e}$ appearing in the bidomain model that can be expressed by

[TABLE]

Exploiting the fact that $D_{e}$ and $D_{i}$ have the same eigenvectors (again defined by the fiber and transmural directions), we can therefore compute the eigenvalues of $K_{0}$ (and analogously of $K_{1}$) as the harmonic mean of the ones of $D_{i}$ and $D_{e}$ as reported, e.g., in [15, Table 8.1]. The values of the main parameters involved in the numerical simulations are reported in Table 1.

The numerical solution of the background problem (2.1) relies on a Newton-Galerkin scheme. The spatial discretization is performed thanks to the P1-finite element space, i.e., the space of the continuous functions over $\Omega$ which are linear polynomials when restricted on each element of the mesh. The temporal discretization is done via an implicit Euler scheme. This leads to solve a nonlinear problem at each timestep, which is performed via a Newton iterative algorithm. More details on the solver can be found in [32], where a thorough convergence analysis is performed.

5.2 Reconstruction of small inclusions

We study the effectiveness of the reconstruction algorithm in identifying the position of small ischemic regions within the cardiac tissue. For, we employ synthetic measurements of the perturbed boundary voltage in presence of ischemias located in different sections of the cardiac tissue: in the left ventricle, in the septum, or in the right ventricle. In each simulation, we consider circular ischemic regions of radius ranging from $1mm$ to $2mm$. In Figure 2 we report the contour plot of the topological gradient $G$ in four different cases, and superimpose a black line representing the boundary of the ischemia that we aim to reconstruct.

From an analysis of Figure 2 we can deduce that the topological gradient always attains its negative minimum value close to the real position of the ischemia, and the accuracy of the localization may vary according to the position of the ischemia. We achieve significant results also in the case of multiple inclusions. When more than one ischemic region is present, the topological gradient shows a local minimum close to each region. As depicted in the captions, some experiments employ the knowledge of the voltage on the epicardium and some on one endocardiac surface. The reconstruction appears fairly accurate in both cases. We remark that the algorithm is always effective when using epicardiac measurements, but fails in detecting ischemic regions located in a ventricle whenever the measurements are acquired on the inner surface of the other one. Nevertheless, even in those cases the technique avoids false positive detection, i.e., the topological gradient does not show significant local minima in wrong locations.

5.3 Reconstruction from noisy measurements

We now focus on the performance of the algorithm in presence of noisy measurements, in order to assess the stability of the algorithm with respect to small perturbations of the boundary data of the form

[TABLE]

where $\eta(x,t)$ is a standard Gaussian random variable for each point $x$ and instant $t$, and $\rho\in[0,1]$ is the noise level.

In Figure 3 we report the contour plot of the topological gradient computed with noiseless measurements and compare it with the ones obtained with growing levels of noise. The algorithm shows to be stable even under significantly corrupted measurements: up to the level $\rho=0.15$, the negative region in correspondence to the exact position is clearly identifiable.

5.4 Reconstruction of large inclusions

We eventually remark that the proposed algorithm produces significant results also when applied to measurements associated to large ischemic regions. The identification of larger regions could be performed by means of iterative reconstruction algorithms, as the one proposed in [7] for a semilinear elliptic problem. Nevertheless, as it is shown in Figure 4, the information coming from the topological gradient could be a suitable initial guess for such iterative algorithms.

6 Conclusions

In this paper, we have studied the identification of small ischemic regions in the cardiac tissue, represented by discontinuous alterations in the conductivity and in the nonlinear reaction term of the monodomain model, taking advantage of the measurement of the voltage on the boundary of the heart. We have extended the existing results regarding the well-posedness of the direct problem and derived a rigorous asymptotic expansion of the perturbed boundary voltage, which allows to formulate an effective reconstruction algorithm based on topological optimization.

We foresee several significant extensions of the presented analysis, both from an analytical and a numerical viewpoint. The coupling of the monodomain model of the heart with a passive conductor model for the surrounding torso (see [11]) would enable us to make use of ECG data for reconstruction purposes, possibly comparing the results with the ones obtained in [23],[38] with a stationary version of the bidomain model.

The effectiveness of the reconstruction algorithm should also be tested in a three-dimensional setting; in such a context, due to the high computational cost of the numerical simulation, a Reduced Order Modeling approach could be considered, as recently analysed in [26].

Finally, the detection of arbitrarily large ischemic regions can be tackled, exploiting an analogous strategy as the one adopted in [6] (possibly reducing the computational effort by employing an adaptive solver of the monodomain system, as proposed in [32]).

Acknowledgments

The authors are grateful to Roberto Gianni for his useful suggestions. The authors thank the New York University of Abu Dhabi for its kind hospitality. The work of the authors has been supported by GNAMPA (Gruppo Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni). The numerical simulations presented in this work have been performed thanks to the MATLAB library redbKIT, [24].