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

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

Elena Beretta Dipartimento di Matematica ”F. Brioschi”, Politecnico di Milano (elena.beretta@polimi.it)    Cecilia Cavaterra Dipartimento di Matematica, Università degli Studi di Milano (cecilia.cavaterra@unimi.it)    M.Cristina Cerutti Dipartimento di Matematica ”F. Brioschi”, Politecnico di Milano” (cristina.cerutti@polimi.it)    Andrea Manzoni CMCS-MATHICSE-SB, Ecole Polytechnique Fédérale de Lausanne (andrea.manzoni@epfl.ch)    Luca Ratti Dipartimento di Matematica ”F. Brioschi”, Politecnico di Milano” (luca.ratti@polimi.it)
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 an inhomogeneity (where the coefficients of the equation are altered) located inside a domain starting from observations of the potential on the boundary . 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.

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 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

(2.1)

where is a bounded set with boundary , and . Here is the domain occupied by the ventricle, is the (transmembrane) electric potential, is a nonlinear term modeling the ionic current flows across the membrane of cardiac cells, is the conductivity tensor of the healthy tissue, and 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 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 , a well-known option would be to replace by where

and solves the following ODE initial value problem, ,

for suitable (constant) parameters , , , . This is the so-called FitzHugh-Nagumo model for the ionic current, and the gating variable 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

(2.2)

where is a parameter determining the rate of change of in the depolarization phase, and are given constant values representing the resting, threshold and peak potentials, respectively. Possible values of the parameters are, e.g., , and , , see [41]. Note that both the sharpness of the wavefront and its propagation speed strongly depend on the value of the parameter .

Consider now a small inhomogeneity located in a measurable bounded domain , such that there exist a compact set , with , and a constant satisfying

(2.3)

Moreover, we assume

(2.4)

In the inhomogeneity the conductivity coefficient and the nonlinearity take different values with respect the ones in . The problem we consider is therefore

(2.5)

where stands for the characteristic function of a set . Here

(2.6)

with .

3 Well posedness of the direct problem

Problem (2.1) thus describes the propagation of the initial activation 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 and we assume that

(3.1)
(3.2)

Moreover, let us set

(3.3)

The following well posedness result holds.

Theorem 3.1.

Assume (2.2), (3.1), (3.2). Then problem (2.1) admits a unique solution such that

(3.4)
(3.5)

where is a positive constant depending (at most) on .

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].

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 such that

(3.6)

Moreover, and the following estimate holds

(3.7)

where is a positive constant depending (at most) on and .

Proof.

Throughout the proof will be as in the statement of the Theorem.

Recalling the definition of , there exist such that

We formulate problem (2.5) in the weak form

(3.8)

Setting , (3.8) becomes

(3.9)

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

(3.10)

the bilinear form is coercive. Indeed

(3.11)

where is a positive constant depending on and .

Through the classical Faedo-Galerkin approximation scheme it is possible to prove that problem (2.5) admits a unique weak solution satisfying (3.6).

In order to obtain further regularity for , let be a sequence such that

and formulate the approximating problems

(3.12)

Using the same arguments as in the proof of Theorem 3.1, we can prove that, , problem (3.12) admits a unique solution such that

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

where are some positive constants independent of .

An application of [38, Thm. 8.1] implies that, up to a subsequence, strongly in , so that a.e. in and . Since problem (2.5) has a unique solution (cf. (3.8)), we conclude that and satisfies

(3.13)

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). ∎

4 Energy estimates for

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

Proposition 4.1.

Assume (2.2), (2.6), (3.1), (3.2). Setting , then

(4.1)
(4.2)

Moreover, there exists such that

(4.3)

Here stands for a positive constant depending (at most) on .

Proof.

Throughout the proof will be as in the statement of the Theorem.

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

(4.4)

where we have set and

(4.5)

being a value between and . By means of (3.4), (3.13) and recalling (3.3), we have

(4.6)

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

Adding and subtracting and applying (3.11) we obtain

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

so that

(4.7)

and finally, see (3.5),

Recalling , an application of Gronwall’s Lemma implies

(4.8)

so that (4.1) follows. Integrating now inequality (4.7) on we get

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

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

(4.9)

Let’s now introduce the auxiliary function , solution to the adjoint problem

(4.10)

By the change of variable , problem (4.10) is equivalent to

(4.11)

where we have set .

Since is bounded in and , standard arguments show that problem (4.11) admits a unique solution such that (see [28, Ch.4])

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

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

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

(4.12)

and then

(4.13)

Recalling , an application of Gronwall’s Lemma gives

(4.14)

Let’s now multiply the first equation in (4.11) by and integrate over . We get

An application of Young’s inequality gives

and then

(4.15)

Combining (4.15) and (4.14), integrating in time on , and using we deduce

so that

(4.16)

The same computations also gives

(4.17)

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

(4.18)

Recalling the definition of and , by estimates (4.16) and (4.18) we get

(4.19)

Finally, we want to prove that there exists such that

(4.20)

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

(4.21)

Moreover, again from (4.19) we have

(4.22)

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

(4.23)

and therefore, using (4.19),

(4.24)

so that (4.20) holds for any .

Let us now multiply the evolution equation in (4.9) by and the evolution equation in (4.10) by , respectively. Integrating on we obtain

(4.25)
(4.26)

Summing up (4.25) and (4.26) we obtain

subsequently, an integration in time on gives

(4.27)

Recalling the conditions at time for and at time for , we get

So that (4.27) becomes

(4.28)

Using now Hölder inequality we deduce

where we may choose for instance and .

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

and therefore

(4.29)

Thanks to (3.5) we also have

Finally, using again Hölder inequality and (4.2), and recalling that , we obtain