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 [art:ECG] and [art:gerbeau]).
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 [art:BCMP],[art:BMR],[art:bccmr]) 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[book:pavarino] and [art:boulakia]) that the monodomain model can be derived from the more complex bidomain one, introduced for the first time in [phd:tung]
, 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.,[monovsbi1],[monovsbi2]).
As a first attempt, we limit ourselves to analyze a coupled system of two equations in the potential and the recovery variable
. 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[book:pavarino] 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 [art:alpanf]), even though the analysis could be extended also to the other two.
In order to detect ischemic regions, we extend the approach of [art:bccmr]
, 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., [book:ammari-kang],[art:cfmv],[art:ammari2012]). A similar approach has also been used in Thermal Imaging (see, e.g., [art:aikk]). Here, we are able to extend in a non trivial way the results obtained previously in [art:BMR] and [art:bccmr] 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 [book:pavarino, Chapter 3] for a general overview. In [art:BCP] 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 [art:veneroni], 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 [pierrecoudiere] 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 [book:smoller] and [book:henry]). 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 [giga].
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 [pierrecoudiere], is proposed in [phd:Luca, 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
where is the region occupied by the hearth tissue,
is the outward unit normal vector to the boundary, and . A slightly different formulation of the monodomain model (see, e.g., [book:pavarino]) 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 . 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 , the perturbed case is described by the model
where is the indicator function of and . We now specify the requirements on the domain, the coefficients and the source terms.
bounded domain, , and ;
the inclusion is well separated from the boundary, i.e.,
, , and on ;
we assume as in the Aliev-Panfilov model, namely:
being and .
The requirements on the matrix-valued functions and are satisfied by the conductivity tensors prescribed in the model under consideration (see [book:pavarino, book:sundes-lines]). 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 the minimum eigenvalue of and by the maximum eigenvalue of .
As an immediate consequence of (2.4), and satisfy the Tangency condition on the rectangle , (see [amann]), i.e., indicating by a generalized outward normal on ( for all , is such that ) then,
Moreover, the functions are Lipschitz continuous on with constants .
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.
Theorem 1 (Unperturbed problem).
The proof of Theorem 1 is derived by means of classical fixed point argument (see [book:pao, Chapter 8, Sections 9 and 11]). A detailed proof is reported in [phd:Luca, 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 [art:BCP] 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).
To start with, note that uniqueness of the weak solution of (2.2) has been shown in the case of Aliev-Panfilov model in [art:kunisch-wagner, Theorem 1.3] as a byproduct of a stability result obtained by exploiting the specific nonlinear expression of and .
We proceed now introducing a sequence of regularized problems 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 is an indicator function, surely ; by density arguments and according to (2.3),
being the space of functions with compact support in . Define and let be the solution of the following problem
where we have used the fact that on . We observe that, for any fixed , 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 and satisfy the Tangency condition on . Hence, from Theorem 1 we deduce , for all .
We now prove that from the convergence of to a weak solution of (2.2) holds. We start by proving some a priori estimates. Consider the weak form of the problem solved by and take the classical solutions , as test functions. Then, we get
Recall now that is the minimum eigenvalue of , whereas is the maximum eigenvalue of . Moreover, since , , are uniformly bounded independently of (indeed, and ) and are continuous, we can introduce and , which are independent of . Hence, by Young inequality,
Using Gronwall’s inequality, we get
Integrating from to (2.9) and using last estimate it also follows that
A bound for the norm of can be found by considering that, for each ,
and computing the norm in time
Analogously, one proves that , with independent of .
As a consequence of the uniform bounds we can ensure that (up to a subsequence) , , , such that
Using the result contained in [book:robinson, Theorem 8.1] and the uniqueness of the weak solution of the perturbed problem we have that (see [book:robinson, Theorem 8.1]), hence a.e. in . Moreover, also a.e. in as it can be straightforwardly obtained by the expression
Consider now the weak form of the problem solved by : , .
Taking the limit in , exploiting the weak convergence of , , , , we obtain that , in and that the limit satisfies, , ,
Indeed, the convergence of the terms involving the time derivatives is a direct consequence of the weak convergence of , 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 and . The convergence of the diffusion term is obtained by combining the weak convergence of in , the pointwise (a.e.) convergence of and the uniform bound on . According to (2.12), we can ensure that the limit is the weak solution of (2.2).
The weak solution is moreover a pointwise (a.e.) limit of the regularized solutions . As a consequence, the uniform bound on is valid also for the limit: a.e. in . This allows to prove the additional Hölder regularity of . Indeed, since , , we can apply Theorem 10.1 of [book:lady, Chapter 3] on the first equation in (2.2)
to get . Now, we can extend the regularity result up to the boundary due to the hypothesis on and on contained in Assumption 1, and conclude . Using the analytic expression of that can be obtained by (2.10) and the regularity of , we can also deduce that . ∎
3 Analysis of the inverse problem
We now tackle the inverse problem of identifying a perturbation from boundary measurements. Suppose to know , the trace of the solution of (2.2) in presence of an unknown inclusion. The inverse problem reads
Although the analysis of the direct problem has been performed in presence of an arbitrary inclusion , 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 [art:BCMP], [art:bccmr]. In particular, we consider a family of inclusions satisfying (2.3) for each and such that
3.1 Energy estimates
We now derive some energy estimates for the difference between and , the solution of the unperturbed problem (2.1) in terms of when .
From now on we will indicate by a positive constant depending on the data, independent of and that may vary also in the same line.
Under Assumption 2.1 the following inequalities hold
, , .
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 , we can deduce
An application of Gronwall’s lemma to (3.4) entails that , and ultimately can be bounded by . This allows to conclude the first two statements.
Observe now that the pair is also the solution of
By the mean value theorem, there exist two pairs of functions , s.t. ***it follows by an application of the Lagrange’s mean value theorem on the real-valued function on the interval ; the same holds for .
By definition, and are convex combinations of and , thus they assume values in the rectangle . Let now be the solution of the adjoint problem
with initial conditions , and homogeneous Neumann boundary condition. Consider the change of variable: , and define , , , and analogously for . Hence and solve
Since and are bounded in , by standard Faedo-Galerkin technique we obtain that the solution of (3.7) exists and is unique with the properties , , , . Moreover
Additional regularity of can be proved with an analogous argument as in [book:lady, Chapter 4, Theorem 9.1], applied on the first equation in (3.7). Indeed, by the regularity of and the square integrability of , we can conclude that , and
Moreover, multiplying the first two equations in (3.7) respectively by and and integrating on , straightforward computations allow to conclude that with