On the inverse problem of detecting cardiac ischemias:
theoretical analysis and numerical reconstruction
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  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.
Mathematical and numerical models of computational electrophysiology can provide quantitative tools to describe electrical heart function and disfunction , 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 .. 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  and  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 . 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  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., ). We use these results to set a reconstruction procedure for detecting the inclusion. To this aim, as in , 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 ).
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
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.,  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.
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 . 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
Moreover, we assume
In the inhomogeneity the conductivity coefficient and the nonlinearity take different values with respect the ones in . The problem we consider is therefore
where stands for the characteristic function of a set . Here
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
Moreover, let us set
The following well posedness result holds.
3.2 Well posedness of the perturbed problem
The well-posedness of the perturbed problem (2.5) is provided by the following theorem.
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
Setting , (3.8) becomes
Observe that, thanks to following the Poincaré type inequality in [7, formula (A.4)]
the bilinear form is coercive. Indeed
where is a positive constant depending on and .
In order to obtain further regularity for , let be a sequence such that
and formulate the approximating problems
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 .
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.
Throughout the proof will be as in the statement of the Theorem.
where we have set and
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
and finally, see (3.5),
Recalling , an application of Gronwall’s Lemma implies
In order to obtain the more refined estimate (4.3), observe that also solves problem
Let’s now introduce the auxiliary function , solution to the adjoint problem
By the change of variable , problem (4.10) is equivalent to
where we have set .
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
Recalling , an application of Gronwall’s Lemma gives
Let’s now multiply the first equation in (4.11) by and integrate over . We get
An application of Young’s inequality gives
The same computations also gives
Finally, we want to prove that there exists such that
To this aim, on account of (4.19) and Sobolev immersion theorems, we deduce
Moreover, again from (4.19) we have
From well-known interpolation estimates (cf. ) we infer
and therefore, using (4.19),
so that (4.20) holds for any .
subsequently, an integration in time on gives
Recalling the conditions at time for and at time for , we get
So that (4.27) becomes
Using now Hölder inequality we deduce
where we may choose for instance and .