On the inverse problem of detecting cardiac ischemias: theoretical analysis and numerical reconstruction
Elena Beretta Cecilia Cavaterra Maria Cristina Cerutti Andrea Manzoni Luca Ratti







equationsection
 16cm 23cm
-15mm -20mm
0pt
theoremTheorem[section]
corollary[theorem]Corollary
definition[theorem]Definition
lemma[theorem]Lemma
proposition[theorem]Proposition
remark[theorem]Remark
#1
    #1(<ref>) Proof. 3mmexdraw
    On the inverse problem of detecting cardiac ischemias: 
     theoretical analysis and numerical reconstruction
        Elena BerettaDipartimento di Matematica "F. Brioschi", Politecnico di Milano ( elena.beretta@polimi.it)  Cecilia CavaterraDipartimento di Matematica, Università degli Studi di Milano ( cecilia.cavaterra@unimi.it) M.Cristina CeruttiDipartimento di Matematica "F. Brioschi", Politecnico di Milano" ( cristina.cerutti@polimi.it) Andrea ManzoniCMCS-MATHICSE-SB, Ecole Polytechnique Fédérale de Lausanne ( andrea.manzoni@epfl.ch) Luca RattiDipartimento di Matematica "F. Brioschi", Politecnico di Milano" ( luca.ratti@polimi.it)
        December 30, 2023 [file: ]
    ===============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
    
    
    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 <cit.> 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. 
    
    
    § INTRODUCTION
    
    
    Mathematical and numerical models of computational electrophysiology can provide quantitative tools to describe electrical heart function and disfunction <cit.>, 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 electrocardiography[The  inverse problem of electrocardiography aims at recovering the epicardial potential (that is, at the heart surface)
    from body surface measurements <cit.>. 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 <cit.>.]. 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 <cit.> and <cit.> 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 <cit.>. Moreover, despite  the inverse problem of ischemia identification from measurements of surface potentials has been tackled in an
    optimization framework for numerical purposes <cit.>, 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 <cit.>.
    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 <cit.> 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.,  <cit.>). A similar approach has also been used in Thermal Imaging (see, e.g., <cit.>).
    
    We use these results to set a reconstruction procedure for detecting the inclusion. To this aim, as in <cit.>,
    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 <cit.>).
    
    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.
    
    
    
    
    
    
    
    
    § 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 <cit.>.
    Throughout the paper we consider the following (background) initial and boundary value problem
    
        ν C_m u_t-  div (k_0∇ u) + ν f(u)  = 0,    in  Ω× (0,T),
           ∂ u/∂ n  = 0,     on   ∂Ω× (0,T), 
            u(0) = u_0,     in  Ω,
        2truecm
    
    where Ω⊂ R^3 is a bounded set with boundary ∂Ω, and
    k_0 ∈ℝ, k_0 >0. Here Ω is the domain occupied by the ventricle, u is the (transmembrane) electric potential, f(u) is a nonlinear term modeling the ionic current flows across the membrane of cardiac cells,
    k_0 is the conductivity tensor of the healthy tissue, C_m >0  and ν >0 are two constant coefficients representing the membrane capacitance and the surface area-to-volume ratio, respectively.
    For the sake of simplicity we deal with an insulated heart, namely we do not consider the effect of the surrounding torso, which behaves as a passive conductor.
    The initial datum u_0 represents the initial activation of the tissue, arising from the propagation of the electrical impulse in the cardiac conduction system. This equation yields a macroscopic model of the cardiac
    tissue, arising from the superposition of intra and extra cellular media, both assumed to occupy the whole heart volume (bidomain model), making the hypothesis that the extracellular and the intracellular conductivities are proportional quantities. Concerning the mathematical analysis of both the monodomain and the bidomain models, some results on the related direct problems have been obtained for instance in <cit.>.
    
    We thus  assume a phenomenological model to describe  the effect of ionic currents through a nonlinear function of the potential.
    We  neglect the coupling with the ODE system modeling the evolution of the so-called  gating variables, which represent the amount of open channels per unit area of the cellular membrane and thus regulate the transmembrane
    currents.
    
    
    
    
    
    
    In the case of a single gating variable w, a well-known option would be to replace  f by g =g(u,w) where
    
        g(u, w) = - β u(u-α)(u-1) - w,
    
    and w solves the following ODE initial value problem, ∀  x ∈Ω,
    
        ∂ w/∂ t = ρ(u- γ w)      (0,T)   ,     w(0) = w_0,
    
    for suitable (constant) parameters β, α, ρ, γ. This is the so-called  FitzHugh-Nagumo model for the ionic current, and the  gating variable w is indeed a  recovery function allowing to take into
    account the depolarization phase. See, e.g., <cit.> for more details. In our case, the model (<ref>) is indeed widely used to characterize the large-scale propagation of the
    front-like solution  in the cardiac excitable medium.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    As suggested in <cit.> and <cit.>, hereon we consider the cubic function
    
    
    
        f(u) = A^2(u-u_1)(u-u_2)(u-u_3),     u_i ∈ℝ,     u_1 < u_2 < u_3,
    
     where   A>0 is a parameter determining the rate of change of u in the depolarization phase, and u_1 <  u_2 < u_3 are given constant values representing the resting, threshold and  peak potentials,
     respectively.  Possible values of the parameters are, e.g., u_1 = -85 mV, u_2 = - 65 mV and u_3 = 40 mV, A=0.04, see <cit.>.
     Note that both the sharpness of the wavefront and its propagation speed  strongly depend on the value of the parameter A.
    
    
    
    
    
    Consider now a small inhomogeneity located in a measurable bounded domain ω_ε⊂Ω, such that there exist a
    compact set K_0, with ω_ε⊂ K_0 ⊂Ω, and a constant d_0 >0 satisfying
    
        dist(ω_ε, Ω\ K_0) ≥ d_0 >0.
    
    Moreover, we assume
    
        |ω_ε| >0,       lim_ε→ 0 |ω_ε| = 0.
    
    In the inhomogeneity ω_ε the conductivity coefficient and the nonlinearity take different values with respect the ones in
    Ω\ω_ε. The problem we consider is therefore
    
        ν C_m u^ε_t-  div (k_ε∇ u^ε) + νχ_Ω\ω_ε f(u^ε) = 0,
            in  Ω× (0,T),
           ∂ u^ε∂ n  = 0,     on   ∂Ω× (0,T), 
            u^ε(0) = u_0,     in  Ω,
        2truecm
    
    where χ_D stands for the characteristic function of a set D ⊂ℝ^3. Here
    
        k_ε = (k_0 - k_1) χ_Ω\ω_ε + k_1 =
           k_0      in  Ω\ω_ε, 
           k_1      in  ω_ε,
    
    with k_0, k_1 ∈ℝ,   k_0 > k_1 >0.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    § WELL POSEDNESS OF THE DIRECT PROBLEM
    
    
    Problem (<ref>)  thus describes the propagation of the initial activation  u_0 in an insulated heart portion (e.g., the left ventricle), and hereon will be referred to as the  background problem;
    we devote Section  <ref> 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 <ref>.
    
    
    
    
     §.§ Well posedness of the background problem
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    For the sake of simplicity, throughout the paper we set
      ν = C_m = 1 and we assume that
     
        Ω∈ C^2+α,  α∈ (0,1),
    
        u_0 ∈ C^2+α(Ω),
             u_1< u_0(x) <u_3   ∀  x ∈Ω,     ∂ u_0 (x)∂ n = 0   ∀  x ∈∂Ω.
    
    
    
    
    Moreover, let us set
    
        M_1:=f_C([u_1,u_3]),      M_2 := f^'_C([u_1,u_3]).
    
    The following well posedness result holds.
    
    Assume (<ref>), (<ref>), (<ref>). Then problem (<ref>) admits a unique solution u∈ C^2+ α, 1 + α/2(Ω× [0,T])
    such that
    
        u_1 ≤ u(x,t) ≤ u_3,    (x,t) ∈Ω× [0,T],
    
        u _C^2+ α, 1 + α/2(Ω× [0,T])≤ C,
    
    where  C is a positive constant depending (at most) on k_0, T, Ω, M_1, M_2,u_0_C^2 + α(Ω).
    
    We omit the details of the proof since (<ref>) can be easily obtained using the results in <cit.>
    and (<ref>) by means of <cit.>.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
     §.§ Well posedness of the perturbed problem
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    The well-posedness of the perturbed problem (<ref>) is provided by the following theorem.
    
    Assume (<ref>), (<ref>), (<ref>), (<ref>). Then problem (<ref>) admits a unique solution u^ε such that
    
        u^ε∈ L^2(0,T; H^1(Ω)) ∩ C([0,T]; L^2(Ω)),    u_t^ε∈ L^2(0,T; (H^1(Ω))^') + L^4/3(Ω× (0,T)).
    
    Moreover, u^ε∈ C^α, α/2(Ω× [0,T]) and the following estimate holds
    
        u^ε_C^α, α/2(Ω× [0,T])≤ C,
    
    where C is a positive constant depending (at most) on k_0, k_1, T, Ω, u_0_C^α(Ω) and M_1.
    
    Throughout the proof C will be as in the statement of the Theorem.
    
    Recalling the definition of f, there exist k ≥ 0,  α_1 >0, α_2 >0,  λ >0 such that
    
        α_1u^4 - k ≤ f(u)u ≤α_2u^4 + k,      f^'(u) ≥ - λ.
    
    We formulate problem (<ref>)  in the weak form
    
        ∫_Ω u^ε_t v + ∫_Ω k_ε∇ u^ε·∇ v
        + ∫_Ωχ_Ω\ω_ε f(u^ε)v = 0 ,   ∀  v ∈ H^1(Ω).
    
    Setting f̃(u) = f(u) - u, (<ref>) becomes
    
        ∫_Ω u^ε_t v + ∫_Ω k_ε∇ u^ε·∇ v
        + ∫_Ωχ_Ω\ω_ε u^ε v
        + ∫_Ωχ_Ω\ω_εf(u^ε)v = 0 ,   ∀  v ∈ H^1(Ω).
    
    Observe that, thanks to following the Poincaré type inequality in <cit.>
        z^2_H^1(Ω)≤ S(Ω)  (∇ z^2_L^2(Ω)
        + z^2_L^2(Ω\ω_ε)),   ∀  z ∈ H^1(Ω),
    
    the bilinear form
    a_ε(u^ε,v) =  ( ∫_Ω k_ε∇ u^ε·∇ v
    + ∫_Ω\ω_ε u^ε v ) is coercive.
    Indeed
    
        a_ε(u^ε,u^ε) =  ∫_Ω k_ε |∇ u^ε|^2 + ∫_Ω\ω_ε (u^ε)^2
        ≥ Su^ε^2_H^1(Ω),
    
    where S is a positive constant depending on Ω and k_1.
    
    Through the classical Faedo-Galerkin approximation scheme it is possible to prove that problem (<ref>) admits a unique weak solution u^ε
    satisfying (<ref>).
    
    In order to obtain further regularity for u^ε, let  {ϕ_n} be a sequence such that
    
        ϕ_n ∈ C^1(Ω),  0 ≤ϕ_n(x) ≤ 1,  ∀  x ∈Ω,  ϕ_n(x) = 1,   ∀  x ∈Ω\ω_ε,  and ϕ_n →χ_Ω∖ω_ε  in  L^∞(Ω),
    
    and formulate the approximating problems
    
        u^n_t-  div (((k_0 - k_1)ϕ_n + k_1)∇ u^n) + ϕ_n f(u^n) = 0,
            in  Ω× (0,T),
           ∂ u^n ∂ n  = 0,     on   ∂Ω× (0,T), 
            u^n(0) = u_0,     in  Ω.
    
    Using the same arguments as in the proof of Theorem <ref>, we can prove that, ∀  n ∈ℕ, problem (<ref>)
    admits a unique solution u^n such that
    
        u^n ∈ C(Ω× [0,T]),    u_1 ≤ u^n(x,t) ≤ u_3,    (x,t) ∈Ω× [0,T].
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    Moreover, by means again of a standard Faedo-Galerkin approximation scheme (for any n) we can prove that the solution to problem (<ref>)
    satisfies
    
        u^n ∈ L^2(0,T; H^1(Ω)) ∩ C([0,T]; L^2(Ω)),
           u_t^n ∈ L^2(0,T; (H^1(Ω))^') + L^4/3(Ω× (0,T)),
    
        u^n^2_L^∞(0,T; L^2(Ω))≤ C,     u^n^2_L^2(0,T; H^1(Ω))≤ C,
    
        u_t^n|^2_L^4/3(0,T; (H^1(Ω))^')≤ C,     ϕ_nf(u^n)^2_L^4/3(Ω× (0,T))≤ C,
    
    where C are some positive constants independent of n.
    
    An application of <cit.> implies that, up to a subsequence,
    u^n →ζ strongly in L^2(Ω× (0,T)),
    so that
    u^n →ζ a.e. in Ω× (0,T) and
    ϕ_nf(u^n) ⇀χ_ω_εf(ζ)   in   L^4/3(Ω× (0,T)).
    
    
    
    
    
    
    Since problem (<ref>) has a unique solution (cf. (<ref>)), we conclude that ζ=u^ε
    and satisfies
    
        u^ε^2_L^∞(0,T; L^2(Ω))≤ C,      u^ε^2_L^2(0,T; H^1(Ω))≤ C,
          u_t^ε^2_L^4/3(0,T; (H^1(Ω))^')≤ C,
           u_1 ≤ u^ε (x,t)≤ u_3,   in  Ω× [0,T].
    
    Considering now the interior regularity result in <cit.> (see also <cit.>)
    and the regularity up to the boundary contained in <cit.>, then we deduce (<ref>).
    
    
    § ENERGY ESTIMATES FOR 𝐮^Ε - 𝐮
    
    
    In this section we prove some energy estimates for the difference between  u^ε and u, solutions to problem (<ref>) and problem (<ref>), respectively, that are crucial to establish the asymptotic formula for u^ε-u of Theorem 4 in Section 5.
    
    Assume (<ref>), (<ref>), (<ref>), (<ref>). Setting w := u^ε - u, then
    
        w_L^∞(0,T;L^2(Ω))≤ C|ω_ε|^1/2,
    
        w_L^2(0,T;H^1(Ω))≤ C|ω_ε|^1/2.
    
    Moreover, there exists 0<β<1 such that
    
        w_L^2(Ω× (0,T))≤ C |ω_ε|^1/2 + β.
    
    Here C stands for a positive constant depending (at most) on k_0, k_1, Ω, T, M_1, M_2, u_0_C^2+α(Ω).
    
    Throughout the proof C will be as in the statement of the Theorem.
    
    On account of the assumptions, Theorems <ref> and <ref> hold. Then w solves the problem
    
        w_t-  div (k_ε∇ w)
        + χ_Ω / ω_ε p_ε w
        = -  div (kχ_ω_ε∇ u) + χ_ω_ε f(u),   in Ω× (0,T),
           ∂ w ∂ n  = 0,   on  ∂Ω× (0,T), 
            w(0) = 0,    in Ω,
    
    where we have set k := k_0 - k_1 > 0 and
    
        p_ε w := f^'(z_ε)w = f(u^ε) - f(u),
    z_ε being a value between u^ε and u. By means of (<ref>), (<ref>)
    and recalling (<ref>), we have
    
        u_1 ≤ z_ε≤ u_3,      |p_ε| = |f^'(z_ε)| ≤ M_2.
    
    Multiplying the first equation by w in (<ref>) by w  and integrating by parts over Ω, we get
    
        1/2d/dt∫_Ω w^2 + ∫_Ω k_ε|∇ w|^2
        + ∫_Ωχ_Ω / ω_ε p_ε w^2
        = ∫_Ωkχ_ω_ε∇ u∇ w + ∫_Ωχ_ω_ε f(u)w.
    
    Adding and subtracting ∫_Ωχ_Ω∖ω_ε w^2 and applying (<ref>) we obtain
    
        1/2d/dt∫_Ω w^2 + Sw^2_H^1(Ω)≤∫_ω_εk∇ u∇ w + ∫_ω_ε f(u)w
        - ∫_Ωχ_Ω / ω_ε (p_ε - 1) w^2.
    
    Recalling (<ref>) and (<ref>), thanks to Young's inequality we deduce
    
        1/2d/dt∫_Ω w^2 + Sw^2_H^1(Ω)≤k (k/2S∫_ω_ε  | ∇ u|^2 + S/2k∫_Ω|∇ w|^2)
        + 1/2∫_ω_ε (f(u))^2 
        + ∫_Ω (M_2 + 3/2) w^2,
    
    so that
    
        1/2d/dt∫_Ω w^2 + S/2w^2_H^1(Ω)≤(k)^2/2S∫_ω_ε  | ∇ u|^2
        + 1/2∫_ω_ε M_1^2
        +  (M_2 + 3/2) ∫_Ω w^2,
    
    and finally, see (<ref>),
    
        d/dtw(t)^2_L^2(Ω)≤ C  (|ω_ε| +  w(t)^2_L^2(Ω) ).
    
    Recalling w(0) = 0, an application of Gronwall's Lemma implies
    
        w(t)^2_L^2(Ω)≤ C|ω_ε|,    t ∈ (0,T),
    
    so that (<ref>) follows.
    Integrating now inequality (<ref>) on (0,T) we get
    
        ∫_Ω w^2(T) + C∫_0^Tw(t)^2_H^1(Ω)dt
        ≤ C  (|ω_ε| +  ∫_0^Tw(t)^2_L^2(Ω)dt),
    
    and a combination with (<ref>) gives (<ref>).
    
    In order to obtain the more refined estimate (<ref>), observe that w also
    solves problem
    
        w_t-  div (k_0∇ w)
        + χ_Ω / ω_ε p_ε w
        = -  div (kχ_ω_ε∇ u^ε) + χ_ω_ε f(u),   in Ω× (0,T),
           ∂ w ∂ n  = 0,   on  ∂Ω× (0,T), 
            w(0) = 0,    in Ω.
    
    Let's now introduce the auxiliary function w, solution to the adjoint problem
    
        w_t +  div (k_0∇w)
        - χ_Ω / ω_ε p_εw
        = - w,   in Ω× (0,T),
           ∂w∂ n  = 0,   on  ∂Ω× (0,T), 
           w(T) = 0,    in Ω.
    
    
    By the change of variable t → T - t, problem (<ref>) is equivalent to
    
        z_t -  div (k_0∇ z)
        + χ_Ω / ω_εp̂_ε z
        = ŵ,   in Ω× (0,T),
           ∂ z ∂ n  = 0,   on  ∂Ω× (0,T), 
            z (0) = 0,    in Ω,
    
    where we have set z(x,t) = w(x, T - t),  p̂_ε (x,t) =  p_ε (x,T - t),  ŵ (x,t)= w(x, T-t).
    
    Since  |χ_Ω / ω_εp̂_ε| is bounded in Ω× (0,T) and w ∈ C^α, α/2(Ω× [0,T]),
    standard arguments show that problem (<ref>) admits a unique solution z such that (see <cit.>)
    
        z ∈ W^2,1_2(Ω× (0,T)) := { z ∈ L^2(Ω× (0,T))   |   z ∈ H^1(0,T; L^2(Ω)) ∩ L^2(0,T;H^2(Ω))}.
    Hereon, up to equation (<ref>), all the equations depend on t and are valid for every t ∈ (0,T); however, we will omit  this dependence for the sake of notation.
    
    Moreover, multiplying the first equation in (<ref>) by z and integrating over Ω, we get
    
        1/2d/dt∫_Ω z^2 +  k_0∫_Ω|∇ z|^2 + k_0∫_Ω z^2
        = ∫_Ωŵ z -∫_Ωχ_Ω / ω_εp̂_ε z^2 + k_0∫_Ω z^2.
    
    By means of Young's inequality and recalling (<ref>), we have
    
        1/2d/dt z_L^2(Ω)^2 +  k_0/2z^2_H^1(Ω)≤1/2k_0ŵ_L^2(Ω)^2  + (M_2 + k_0)  z_L^2(Ω)^2,
    
    and then
    
        d/dt z_L^2(Ω)^2
        ≤1/k_0ŵ_L^2(Ω)^2 + 2(M_2 + k_0)  z_L^2(Ω)^2.
    
    Recalling z(x,0)= 0, an application of Gronwall's Lemma gives
    
        z _L^2(Ω)^2
        ≤ Cŵ_L^2(Ω)^2.
    
    Let's now multiply the first equation in (<ref>) by z_t and integrate over Ω. We get
    
        ∫_Ω z_t^2 +  k_0/2d/dt∫_Ω |∇ z|^2
        = ∫_Ωŵ z_t -∫_Ωχ_Ω / ω_εp̂_ε z z_t.
    
    An application of Young's inequality gives
    
        1/2∫_Ω z_t^2 +  k_0/2d/dt∫_Ω|∇ z|^2
        ≤∫_Ω (ŵ)^2 +∫_Ωχ_Ω / ω_ε (p̂_ε)^2 z^2,
    
    and then
    
        1/2z_t_L^2(Ω)^2 +  k_0/2d/dt∇ z_L^2(Ω)^2
        ≤ŵ_L^2(Ω)^2 +M_2^2 z_L^2(Ω)^2.
    
    Combining (<ref>) and (<ref>),
    
    
    
    integrating in time on (0,t), and using ∇ z(0) = 0 we deduce
    
        ∇ z(t)_L^2(Ω)^2 ≤ C ŵ_L^2(Ω× (0,t))^2,    t ∈ (0,T),
    
    so that
    
        z_L^∞(0,T;H^1(Ω))^2 ≤ C ŵ_L^2(Ω× (0,T))^2.
    
    The same computations also gives
    
    
        z_t_L^2(Ω× (0,T))^2
        ≤ Cŵ_L^2(Ω× (0,T))^2.
    
    Then, an application of standard elliptic regularity results to problem (<ref>) implies (see <cit.>)
    
        z_L^2(0,T;H^2(Ω)^2
        ≤ Cŵ_L^2(Ω× (0,T))^2.
    
    
    Recalling the definition of z and ŵ, by estimates (<ref>) and (<ref>)
    we get
    
        w^2_L^∞(0,T;H^1(Ω)) + w_L^2(0,T;H^2(Ω))^2
        ≤ Cw_L^2(Ω× (0,T))^2,
    
    Finally, we want to prove that there exists p >2 such that
    
        w_L^p(Ω× (0,T)) + ∇w_L^p(Ω× (0,T))≤ Cw_L^2(Ω× (0,T)).
    
    To this aim, on account of (<ref>) and Sobolev immersion theorems, we deduce
    
        w^2_L^6(Ω× (0,T))≤ Cw^2_L^∞(0,T;H^1(Ω))≤ Cw_L^2(Ω× (0,T))^2.
    
    Moreover, again from (<ref>) we have
    
        ∇w∈ L^∞(0,T; L^2(Ω)) ∩ L^2(0,T;L^6(Ω)).
    
    From well-known interpolation estimates (cf. <cit.>) we infer
    
        ∇w^10/3_L^10/3(Ω× (0,T))≤ C∇w^2_L^2(0,T;L^6(Ω))∇w^4/3_L^4/3(0,T;L^2(Ω))
    
    and therefore, using (<ref>),
    
        ∇w^10/3_L^10/3(Ω× (0,T))≤ C w ^2_L^2(Ω× (0,T)) w ^4/3_L^2(Ω× (0,T))≤ C  w ^10/3_L^2(Ω× (0,T)),
    
    so that (<ref>) holds for any p ∈ (2, 10/3].
    
    Let us now multiply the evolution equation in (<ref>) by w and the evolution equation in (<ref>) by w, respectively.
    Integrating on Ω we obtain
    
        ∫_Ω w_t w + k_0 ∫_Ω∇ w ·∇w
        + ∫_Ωχ_Ω / ω_ε p_ε w w
        =  k∫_ω_ε∇ u^ε∇w  + ∫ _ω_ε f(u) w,
    
        ∫_Ωw_t  w -  k_0 ∫_Ω∇w·∇  w
        - ∫_Ωχ_Ω / ω_ε p_εw  w
        = -  ∫ _Ω w^2.
    
    Summing up (<ref>) and (<ref>) we obtain 
    
        ∫_Ω (w_t w + w_t  w) =
        k∫_ω_ε∇ u^ε∇w  + ∫ _ω_ε f(u) w
         -  ∫ _Ω w^2,
    
    subsequently, an integration in time on (0,T) gives
    
        ∫_0^T ∫ _Ω w^2
        = -∫_0^T∫_Ω (w_t w + w_t  w)
        + k∫_0^T∫_ω_ε∇ u^ε∇w  + ∫_0^T∫ _ω_ε f(u) w.
    
    Recalling the conditions at time t=0 for w and at time t=T for w, we get
    
        [  ∫_0^T∫_Ω (w_t w + w_t  w) = ∫_Ω∫_0^T  (w_t w +w_t  w )                                                        ;                                                         = ∫_Ω ((ww)(T) - (ww)(0) - ∫_0^T (w w_t  + w_t  w) )
        =0 ]
    
    So that (<ref>) becomes
    
        ∫_0^T ∫ _Ω w^2
        = k∫_0^T∫_ω_ε∇ u^ε∇w  + ∫_0^T∫ _ω_ε f(u) w.
    
    Using now Hölder inequality we deduce
    
        w_L^2(Ω× (0,T))^2
        ≤∇ u^ε_L^q(ω_ε× (0,T))∇w_L^p(ω_ε× (0,T))
        +f(u)_L^q(ω_ε× (0,T))w_L^p(ω_ε× (0,T)),
    
    where we may choose for instance p=10/3 and q = 10/7.
    
    By means of (<ref>) and  (<ref>),  from the previous inequality we get
    
        w^2_L^2(Ω× (0,T))≤ C w_L^2(Ω× (0,T)) (∇ u^ε_L^q(ω_ε× (0,T))
        + |ω_ε|^1/q),
    
    and therefore
    
        w_L^2(Ω× (0,T))≤ C (∇ u^ε_L^q(ω_ε× (0,T))
        + |ω_ε|^1/q ).
    
    Thanks to (<ref>) we also have
    
        [                              ∇ u^ε_L^q(ω_ε× (0,T)) ≤∇ u^ε - ∇ u_L^q(ω_ε× (0,T)) + ∇ u_L^q(ω_ε× (0,T));                                                                   ≤∇ w _L^q(ω_ε× (0,T)) + C|ω_ε|^1/q. ]
    
    Finally, using again Hölder inequality and (<ref>), and recalling that q ∈ [10/7,)2, we obtain
    
        ∇ w _L^q(ω_ε× (0,T))   ≤ ( ∫_0^T (∫_ω_ε|∇ w (t)|^q 2/p^')^q/2 (∫_ω_ε 1 )^2-q/2 )^1/q
           ≤ |ω_ε|^1/q - 1/2 (∫_0^T∇ w(t)^q_L^2(Ω) )^1/q≤  |ω_ε|^1/q - 1/2∇ w_L^q(0,T; L^2Ω))
           ≤ C(Ω) |ω_ε|^1/q - 1/2∇ w_L^2(Ω× (0,T))≤ C |ω_ε|^1/q.
    
    Combining the previous estimate with (<ref>), since 1/q∈(1/2,7/10] we can conclude that (<ref>) holds with β∈(0,1/5].
    
    
    § THE ASYMPTOTIC FORMULA
    
    
    In this section we derive and prove an asymptotic representation formula for w = u_ε - u in analogy with <cit.> and <cit.>.
    Let Φ =Φ(x,t) be any solution of
    
        Φ_t + k_0ΔΦ  - f^'(u)Φ  = 0,    in  Ω× (0,T),
           Φ(T) = 0,     in  Ω.
        2truecm
    
    
    Our main result is the following
    
    Assume (<ref>), (<ref>), (<ref>), (<ref>). Let u^ε and u be the solutions to (<ref>) and (<ref>)
    and Φ a solution to (<ref>), respectively. Then, there exist a sequence ω_ε_n satisfying (<ref>) and (<ref>)
    with |ω_ε_n|→ 0, a regular Borel measure μ and a symmetric matrix ℳ with elements
    ℳ_ij∈ L^2(Ω, dμ) such that, for ε→ 0,
    
        ∫_0^T∫_∂Ω k_0∂Φ/∂ n(u^ε - u)
        = |ω_ε_n| {∫_0^T∫_Ωk̃ℳ∇ u·∇Φ dμ
        +  ∫_0^T∫_Ω  f(u)Φ dμ + o(1) }.
    
    To prove Theorem <ref>, we need to state some preliminary results. Let  v_ε^(j) and v^(j) be the variational solutions
    (depending only on x ∈Ω) to the problems
    
        (PV_ε)
        {[                                  div (k_ε∇ v_ε^(j))  = 0,   in Ω, ∂ v_ε^(j)/∂ n  = n_j,    on  ∂Ω,                ∫_∂Ω v_ε^(j) = 0, ].
          
        (PV_0) {[                                div (k_0∇ v^(j))  = 0,   in Ω, ∂ v^(j)/∂ n  = n_j,    on  ∂Ω,                 ∫_∂Ω v^(j) =0, ].
    n_j being the j-th coordinate of the outward normal to ∂Ω. It can be easily verified that
    
        v^(j) = x_j - 1/|∂Ω|∫_∂Ωx_j.
    
    The following results hold
    
    Let v_ε^(j) and v^(j) solutions to (<ref>), then there exists C(Ω)>0 such that
    
        v_ε^(j) - v^(j)_H^1(Ω)≤ C(Ω)|ω_ε|^1/2.
    
    Moreover, for some η∈ (0, 1/2), there exists C(Ω, η)>0 such that
    
        v_ε^(j) - v^(j)_L^2(Ω)≤ C(Ω, η)|ω_ε|^1/2 + η.
    
    See Lemma 1 in <cit.>.
    
    Let u and u_ε be the solutions to problems (<ref>) and (<ref>), respectively.
    Consider v_ε^(j) and  v^(j) as in (<ref>).
    Then, for any Φ∈ C^1(Ω× [0,T]) with Φ(x,T) = 0, the folllowing holds as ε→ 0,
    
        ∫_0^T∫_Ω1/|ω_ε|χ_ω_ε∇ u·∇ v_ε^(j)Φ dxdt =
        ∫_0^T∫_Ω1/|ω_ε|χ_ω_ε∇ u_ε·∇ v^(j)Φ dxdt + o(1).
    
    We follow the ideas in <cit.> and <cit.>.
    Since w = u_ε - u, then we obtain the identity
    
        ∫_Ω k_0∇ w·∇ (v^(j)Φ)    = ∫_Ω k_0∇ w·∇ v^(j)Φ
        + ∫_Ω k_0∇ w·∇Φ v^(j)
           = - ∫_Ω k_0 w ∇ v^(j)·∇Φ + ∫_∂Ω k_0 w n_jΦ
        +∫_Ω k_0∇ w·∇Φ v^(j).
    
    Moreover, we have
    
        ∫_0^T∫_Ω k_ε∇ w·∇ (v_ε^(j)Φ)
        = ∫_0^T∫_Ω (  k_ε∇ w·∇ v_ε^(j)Φ
        +  k_ε∇ w·∇Φ v^(j)
        +  k_ε∇ w·∇Φ (v^(j)_ε - v^(j)) )
    
        = ∫_0^T ∫_Ω  k_ε∇ w·∇ v_ε^(j)Φ
        + ∫_Ω k_0∇ w·∇Φ v^(j)
        +  ∫_Ω  (k_ε - k_0)∇ w·∇Φ v^(j) + ∫_Ω k_ε∇ w·∇Φ (v^(j)_ε - v^(j))
    
        = ∫_0^T ( - ∫_Ω k_ε w ∇ v_ε^(j)·∇Φ
        + ∫_∂Ωk_0 w n_jΦ
        + ∫_Ω k_0∇ w·∇Φ v^(j)
    
        + ∫_Ω (k_ε - k_0)∇ w·∇Φ v^(j)
        + ∫_Ω k_ε∇ w·∇Φ (v^(j)_ε - v^(j)) )
    
        = ∫_0^T ( - ∫_Ω k_ε w ∇ v^(j)·∇Φ
        + ∫_∂Ωk_0 w n_jΦ
        + ∫_Ω k_0∇ w·∇Φ v^(j)
    
        + ∫_ω_ε (k_1 - k_0)∇ w·∇Φ v^(j) + ∫_Ω k_ε∇ w·∇Φ (v^(j)_ε - v^(j))
        - ∫_Ω k_ε w ∇ (v^(j)_ε - v^(j))·∇Φ).
    
    A combination with (<ref>) gives
    
        ∫_0^T∫_Ω k_ε∇ w·∇ (v_ε^(j)Φ )
        = ∫_0^T ( ∫_Ω k_0∇ w·∇  (v^(j)Φ )  +  ∫_Ω(k_0 - k_ε) w ∇ v^(j)·∇Φ
    
        + ∫_ω_ε (k_1 - k_0)∇ w·∇Φ v^(j) + ∫_Ω k_ε∇ w·∇Φ (v^(j)_ε - v^(j))
        - ∫_Ω k_ε w·∇ (v^(j)_ε - v^(j))·∇Φ)
    
        = ∫_0^T ( ∫_Ω k_0∇ w·∇ (v^(j)Φ) + ∫_ω_ε (k_1 - k_0)∇ w·∇Φ v^(j)
    
        +  ∫_ω_ε(k_0 - k_1) w ∇ v^(j)·∇Φ + ∫_Ω k_ε∇ w·∇Φ (v^(j)_ε - v^(j))
        - ∫_Ω k_ε w ∇ (v^(j)_ε - v^(j))·∇Φ).
    
    Then, on account of (<ref>), (<ref>), (<ref>), (<ref>) and Schwarz inequality, we get
    
    
    
    
    
    
    
    
    
        ∫_0^T∫_Ω k_ε∇ w·∇  (v_ε^(j)Φ )
        = ∫_0^T ( ∫_Ω k_0∇ w·∇ (v^(j)Φ ) - ∫_ω_εk∇ w·∇Φ v^(j) )
         +o(|ω_ε|).
    
    Let us consider now problem (<ref>). Multiplying both sides of the first equation by v_ε^(j)Φ, and integrating by parts
    on Ω× (0,T), we obtain
    
        ∫_0^T∫_Ω w_tv_ε^(j)Φ + ∫_0^T∫_Ω k_ε∇ w ·∇ (v_ε^(j)Φ)
        + ∫_0^T∫_Ωχ_Ω / ω_ε p_ε w v_ε^(j)Φ
        
        = ∫_0^T∫_ω_εk∇ u·∇(v_ε^(j)Φ) + ∫_0^T∫_ω_ε f(u)v_ε^(j)Φ.
    
    On the other hand, multiplying the first equation in (<ref>) by v^(j)Φ and integrating by parts on Ω× (0,T), we get
    
        ∫_0^T∫_Ω w_t v^(j)Φ + ∫_0^T∫_Ω k_0∇ w ·∇(v^(j)Φ)
        + ∫_0^T∫_Ωχ_Ω / ω_ε p_ε w v^(j)Φ
        
        =  ∫_0^T∫_ω_εk∇ u^ε·∇(v^(j)Φ) + ∫_0^T∫_ω_ε f(u) v^(j)Φ.
    
    A combination of (<ref>), (<ref>) and (<ref>) gives, for ε→ 0,
    
        ∫_0^T∫_ω_εk∇ u·∇(v_ε^(j)Φ) + ∫_0^T∫_ω_ε f(u)v_ε^(j)Φ
        - ∫_0^T∫_Ω w_tv_ε^(j)Φ
        - ∫_0^T∫_Ωχ_Ω / ω_ε p_ε w v_ε^(j)Φ
    
        =  ∫_0^T∫_ω_εk∇ u^ε·∇(v^(j)Φ) + ∫_0^T∫_ω_ε f(u) v^(j)Φ
        - ∫_0^T∫_Ω w_t v^(j)Φ
        - ∫_0^T∫_Ωχ_Ω / ω_ε p_ε w v^(j)Φ
    
        - ∫_0^T∫_ω_εk̃∇ w·∇Φ v^(j)
         +o(|ω_ε|),
    
    from which we deduce
    
        ∫_0^T∫_ω_εk∇ u·∇(v_ε^(j)Φ)
        =  ∫_0^T∫_Ω w_t(v_ε^(j) - v^(j))Φ
        + k∫_0^T∫_ω_ε ( ∇ u^ε·∇(v^(j)Φ)
        - ∇ u^ε·∇Φ v^(j) )
    
        + ∫_0^T∫_ω_εk∇ u·∇Φ v^(j)
        - ∫_0^T∫_Ωχ_Ω / ω_ε p_ε w (v_ε^(j) - v^(j))Φ
        + ∫_0^T∫_ω_ε f(u) (v^(j) - v_ε^(j))Φ
        + o(|ω_ε|).
    
    By means of  (<ref>), (<ref>), (<ref>) and (<ref>), and recalling also (<ref>) and (<ref>),
    an application of the Hölder inequality both in space and time gives, for ε→ 0,
    
        k∫_0^T∫_ω_ε∇ u·∇(v_ε^(j)Φ)
        =  ∫_0^T∫_Ω w_t(v_ε^(j) - v^(j))Φ
        + k∫_0^T∫_ω_ε ( ∇ u^ε·∇ v^(j)Φ
        + ∇ u ·∇Φ v^(j) )
         +o(|ω_ε|)
    
    and then
    
        k∫_0^T∫_ω_ε∇ u·∇ v_ε^(j)Φ
        =  ∫_0^T∫_Ω w_t(v_ε^(j) - v^(j))Φ
        
        + k∫_0^T∫_ω_ε∇ u^ε·∇ v^(j)Φ
        + k∫_0^T∫_ω_ε∇ u ·∇Φ (v^(j)- v_ε^(j))
        +o(|ω_ε|) 
        
        =  ∫_0^T∫_Ω w_t(v_ε^(j) - v^(j))Φ
        + k∫_0^T∫_ω_ε∇ u^ε·∇ v^(j)Φ
        +o(|ω_ε|).
    
    Consider the first term in the last line of (<ref>). Integrating by parts in time and recalling that
    Φ(T) = 0, w(0) = 0, (v^(j)- v_ε^(j))_t = 0, we finally have
    (cf. also (<ref>), (<ref>)), for ε→ 0,
    
        ∫_0^T∫_Ω w_t(v_ε^(j) - v^(j))Φ =
        ∫_Ω [w(v_ε^(j) - v^(j))Φ](T)
        - ∫_Ω [w(v_ε^(j) - v^(j))Φ](0) 
           -  ∫_0^T∫_Ω w(v_ε^(j) - v^(j))_tΦ  -  ∫_0^T∫_Ω w(v_ε^(j) - v^(j))Φ _t
        =  -  ∫_0^T∫_Ω w(v_ε^(j) - v^(j))Φ _t = o(|ω_ε|).
    
    Combining (<ref>) and (<ref>) we get
    
        k∫_0^T∫_ω_ε∇ u·∇ v_ε^(j)Φ
        =  k∫_0^T∫_ω_ε∇ u^ε·∇ v^(j)Φ
        +o(|ω_ε|),   ε→ 0,
    
    then formula (<ref>) is true.
     Proof of Theorem <ref>.
    Following <cit.>, there exist a regular Borel measure μ, a symmetric matrix ℳ with elements
    ℳ_ij∈ L^2(Ω, dμ), a sequence ω_ε_n with |ω_ε_n| → 0 such that
    
        |ω_ε_n|^-1χ_ω_ε_n dx → dμ,
             |ω_ε_n|^-1χ_ω_ε_n∂ v_ε_n^(j)/∂ x_i dx →ℳ_ij dμ,
    
    in the weak^* topology of C^0(Ω).
    On account of  (<ref>), we deduce also
    
        |ω_ε_n|^-1χ_ω_ε_n∂ u(t)/∂ x_i∂ v_ε_n^(j)/∂ x_i dx
        →ℳ_ij∂ u(t)/∂ x_i dμ,   ∀  t ∈ (0,T),
    
    in the weak^* topology of C^0(Ω).
    Moreover, recalling (<ref>), (<ref>) and (<ref>), we get
    
        | ∫_0^T∫_Ωχ_ω_ε_n/|ω_ε_n|∂ u^ε_n/∂ x_i∂ v^(j)/∂ x_i |
        ≤ C,
    
    
    
    
    
    
    
    
    
    
    where C is independent of ε_n.
    Hence
    
        |ω_ε_n|^-1χ_ω_ε_n∂ u^ε_n/∂ x_i∂ v^(j)/∂ x_i dxdt → dν_j
    
    in the weak^* topology of C^0(Ω× [0,T]).
    Combining (<ref>), (<ref>) and (<ref>) we obtain
    
        dν_j = ℳ_ij∂ u(t)/∂ x_i dμ,   ∀  t ∈ (0,T).
    
    
    Now multiply the first equation in (<ref>) by w and the first equation in (<ref>) by Φ on Ω× (0,T).
    
    
    
    Then, integrating by parts, we get
    
        ∫_0^T∫_ΩΦ_tw + ∫_0^T∫_Ω k_0∇Φ·∇ w  - ∫_0^T∫_Ω f^'(u)Φ w
        +∫_0^T∫_∂Ω k_0∂Φ/∂ nw= 0,
    
        ∫_0^T∫_Ω w_tΦ  + ∫_0^T∫_Ω k_0∇ w ·∇Φ + ∫_0^T∫_Ωχ_Ω / ω_ε p_ε w Φ
        =   ∫_0^T∫_ω_εk∇ u^ε·∇Φ +  ∫_0^T∫_ω_ε f(u)Φ.
    
    Summing up the two previous equations, we have
    
        ∫_0^T∫_Ω (w_tΦ + Φ_tw) - ∫_0^T∫_Ω f^'(u)Φ w
        +∫_0^T∫_∂Ω k_0∂Φ/∂ nw + ∫_0^T∫_Ωχ_Ω / ω_ε p_ε w Φ
        
        =   ∫_0^T∫_ω_εk∇ u^ε·∇Φ +  ∫_0^T∫_ω_ε f(u)Φ.
    
    Observe that the following identities hold
    
        ∫_0^T∫_Ω (w_tΦ + Φ_tw) = ∫_Ω ( Φ(T)w(T) - Φ(0)w(0) ) - ∫_0^T∫_ΩΦ w_t +   ∫_0^T∫_ΩΦ w_t =0,
    
    and then, from (<ref>) we infer
    
        ∫_0^T∫_Ω ( χ_Ω / ω_ε p_ε w Φ - f^'(u)Φ w )
        +∫_0^T∫_∂Ω k_0∂Φ/∂ nw
        =   ∫_0^T∫_ω_εk∇ u^ε·∇Φ +  ∫_0^T∫_ω_ε f(u)Φ.
    
    Moreover, on account of (<ref>), we have
    
        ∫_0^T∫_Ω ( χ_Ω / ω_ε p_ε w Φ - f^'(u)Φ w  )
           =∫_0^T∫_Ω ( χ_Ω / ω_ε p_ε w Φ - χ_Ω / ω_ε f^'(u)Φ w  )
        - ∫_0^T∫_ω_ε f^'(u)Φ w
           = ∫_0^T∫_Ωχ_Ω / ω_ε (p_ε -  f^'(u)) w Φ + o(|ω_ε|) = o(|ω_ε|).
    
    The last equality in (<ref>) is a consequence of the regularity of f (see (<ref>), from which |p_ε - f^'(u)| ≤ C|w| follows)
    and (<ref>).
    
    
    
    
    
    
    
    
    
    
    Combining (<ref>) and (<ref>)
    we obtain
    
        ∫_0^T∫_∂Ω k_0∂Φ/∂ nw
        
        = |ω_ε|∫_0^T∫_Ωk |ω_ε|^-1χ_ω_ε∇ u^ε·∇Φ
        + |ω_ε| ∫_0^T∫_Ωχ_ω_ε|ω_ε|^-1 f(u)Φ
        + o(|ω_ε_n|).
    
    And finally, by means of (<ref>), (<ref>) and (<ref>), the formula (<ref>) holds. 
    We would like to emphasize that, with minor changes, the asymptotic expansion extends  to  the case of piecewise smooth anisotropic conductivities
    of the form
    
        𝕂_ε = {[   𝕂_0 Ω∖ω_ε;   𝕂_1   ω_ε ].
    
    where 𝕂_0, 𝕂_1∈ C^∞(Ω) are symmetric matrix valued functions satisfying
    
        α_0|ξ|^2≤ξ^T𝕂_0(x)ξ≤β_0 |ξ|^2,   α_1|ξ|^2≤ξ^T𝕂_1(x)ξ≤β_1|ξ|^2,   ∀ ξ∈ℝ^3, ∀  x ∈Ω,
    
    with 0<α_1<β_1< α_0 < β_0. Then, the asymptotic formula (<ref>) becomes
    
        ∫_0^T∫_∂Ω𝕂_0∇Φ· n (u^ε - u) =
        |ω_ε|∫_0^T∫_Ω( ℳ_i j(𝕂_0-𝕂_1)_ik∂ u/∂ x_k∂Φ/∂ x_j
        +f(u)Φ )dμ  +o(|ω_ε|)
    
    where Φ solves
    
        Φ_t +  div (𝕂_0 ∇Φ) - f^'(u)Φ  = 0,    in  Ω× (0,T),
           Φ(T) = 0,     in  Ω,
        			2truecm
    
    and u is the background solution of
    
        u_t -  div (𝕂_0 ∇ u) + f(u)  = 0,    in  Ω× (0,T),
           𝕂_0∇ u · n = 0,     on  ∂Ω× (0,T), 
            u(0) = 0,     in  Ω.
        			2truecm
    
    The matrix ℳ is called  the polarization tensor associated to the inhomogeneity ω_ε.
    Indeed, all the results of the previous sections can be extended to the case of constant anisotropic coefficients using for instance
    the regularity results contained in <cit.>.
    
    
    § A RECONSTRUCTION ALGORITHM
    
    
    We now use the asymptotic expansion derived in the previous section to set a reconstruction procedure
    for the inverse problem of detecting a spherical inhomogeneity ω_ε from boundary measurements of the potential.
    Following the approach of <cit.>,
    <cit.>, but taking into account the time-dependence of the problem, we introduce the mismatch functional
    
        J(ω_ε) = 1/2∫_0^T ∫_∂Ω(u^ε - u_meas)^2 ,
    
    being u^ε the solution of the perturbed problem (<ref>) in presence of an inclusion ω_ε satisfying hypotheses
    (<ref>), (<ref>).
    It is possible to reformulate the inverse problem in terms of the following minimization problem
    
        J(ω_ε) →min
    
    among all the small inclusions, well separated from the boundary. We introduce the following additional assumption on the exact inclusion
    
        ω_ε = z + ε B = {x ∈Ω s.t.  x = z+ε b, b ∈ B},
    
    being z ∈Ω and B an open, bounded, regular set containing the origin. We remark that we prescribe the geometry of the inclusion to be fixed
    throughout the whole observation time.
    The restriction of the functional J to the class of inclusions satisfying (<ref>) is denoted by j(ε;z).
    We can now define the topological gradientG: Ω→ as the first order term appearing in the asymptotic expansion of the cost
    functional with respect to ε, namely
    
        j(ε;z) = j(0) + |ω_ε| G(z) + o(|ω_ε|),   ε→ 0,
    
    where j(0) = ∫_0^T ∫_∂Ω(u - u_meas)^2 and u is the solution of the unperturbed problem (<ref>).
    Under the assumption that the exact inclusion has small size and satisfies hypothesis
    (<ref>), a reconstruction procedure consists in identifying the point z̅∈Ω where the topological gradient G attaints its minimum.
    Indeed, the cost functional achieves the smallest value when it is evaluated in the center of the exact inclusion. Thanks to the hypothesis of small size, we expect the reduction of the cost functional j to be correctly described by the first order term G, up to a reminder which is negligible with respect to ε.
    
    Nevertheless, in order to define a reconstruction algorithm, we need to efficiently evaluate the topological gradient G. According to the definition,
    
        G(z) = lim_ε→ 0j(ε;z)-j(0)/|ω_ε|.
    
    Evaluating G in a single point z∈Ω would require to solve the direct problem several times in presence of inclusions centered at z
    with decreasing volume.
    This procedure can be indeed avoided thanks to a useful representation formula that can be deduced from the asymptotic expansion (<ref>).
    To show this we need the following preliminary Proposition the proof of which is given in the Appendix.
    
    	Consider the problem
    	
        Φ_t + k_0ΔΦ - f^'(u)Φ  = 0,    in  Ω× (0,T),
           ∂Φ/∂ n  = u^ε - u,     on   ∂Ω× (0,T), 
           Φ(T) = 0,     in  Ω.
        			2truecm
    
    	Given a compact set K ⊂Ω such that d(K,∂Ω) ≥ d_0 > 0 the following estimate holds
    	
        Φ_L^1(0,T;W^1,∞(K))≤ C u^ε - u _L^2(0,T;L^2(∂Ω)).
    
    On account of Proposition <ref>, we deduce the following representation of the topological gradient
    
    The topological gradient of the cost functional j(ε, z) can be expressed by
    
        G(z) = ∫_0^T  ( kℳ∇ u(z) ·∇ W(z)  +  f(u(z))W(z) ),
    
    where W is the solution of the adjoint problem:
    
        W_t + k_0Δ W  - f^'(u)W  = 0,    in  Ω× (0,T),
           ∂ W/∂ n = u - u_meas,     on  ∂Ω× (0,T), 
            W(T) = 0,     in  Ω.
    
    
    Consider the difference
    
        j(ε;z) - j(0)
           = 1/2u^ε - u_meas_L^2(0,T;L^2(∂Ω))^2 - 1/2u - u_meas_L^2(0,T;L^2(∂Ω))^2 
            =  ∫_0^T ∫_∂Ω (u^ε - u)(u - u_meas)dt  + 1/2u^ε - u _L^2(0,T;L^2(∂Ω))^2.
    
    According to (<ref>) and to the definition of the adjoint problem (<ref>), we can express
    
        ∫_0^T ∫_∂Ω (u^ε - u)(u - u_meas)dt =
        |ω_ε| {∫_0^T∫_Ωkℳ∇ u·∇ W dμ
        +  ∫_0^T∫_Ω  f(u)W dμ + o(1) }.
    
    Since we assume (<ref>), the measure μ associated to the inclusion is the Dirac mass δ_z centered in point z (see <cit.>).
    Hence
    
        ∫_0^T ∫_∂Ω (u^ε - u)(u - u_meas)dt =    
        |ω_ε| {∫_0^T kℳ∇ u(z)·∇ W(z)  . 
            + . ∫_0^T f(u(z))W(z) } + o(|ω_ε|).
    
    Moreover, by (<ref>), the second term in the left-hand side of (<ref>) can be expressed as
    
        ∫_0^T ∫_∂Ω (u^ε - u)(u^ε - u)dt =
        |ω_ε| {∫_0^Tkℳ∇ u(z)·∇Φ(z) +  ∫_0^T f(u(z))Φ(z) } + o(|ω_ε|),
    
    where Φ is the solution to (<ref>). Thanks to regularity results on u (see Theorem <ref>) and using Proposition <ref> with
    K = Ω_d_0 = {x ∈Ω s.t. d(x,∂Ω) ≥ d_0}, we obtain
    
        ∫_0^T ∫_∂Ω (u^ε - u)(u^ε - u)dt ≤ C |ω_ε| {∫_0^T | ∇Φ(z)|
        +  ∫_0^T |Φ(z)| } + o(|ω_ε|) 
           ≤ C|ω_ε| u^ε - u _L^2(0,T,L^2(∂Ω)) + o(|ω_ε|)
         ≤ C|ω_ε| u^ε - u _L^2(0,T,H^1(Ω)) + o(|ω_ε|) 
           ≤ C|ω_ε|^3/2 + o(|ω_ε|) = o(|ω_ε|).
    
    Replacing (<ref>) and (<ref>) in (<ref>), we finally get
    
        j(ε;z) - j(0)
        = |ω_ε| {∫_0^T kℳ∇ u(z)·∇ W(z) +  ∫_0^T f(u(z))W(z) } + o(|ω_ε|).
    
    
    
    
    
    
    
    
    
    
    
    
    
    Thanks to the representation formula (<ref>), evaluating the topological gradient of the cost functional requires just the solution
    of two initial and boundary value problems. This yields the definition of
    a one-shot algorithm for the identification of the center of a small inclusion satisfying hypotesis (<ref>) (see Algorithm <ref>).
    
    
    
    Inspired by the electrophysiological application, we consider moreover the possibility to have  partial boundary measurements, i.e.
    the case where the support of the function u_meas is not the whole
    boundary ∂Ω but only a subset Γ⊂∂Ω. In this case, it is possible to formulate a slightly different optimization problem,
    in which we aim at minimizing the mismatch between the measured and the perturbed data just on the portion Γ of the boundary. The same reconstruction algorithm can be devised
    for this problem, by simply changing the expression of the Neumann condition of the adjoint problem.
    
    
    
    
    
    
    	
    
    
    
    	
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    	
    
    
    
    
    
    
    
    
    	
    	
    		
    			
    			
    		
    	
    	
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    § NUMERICAL RESULTS
    
    
    In order to implement Algorithm 1 for the detection of inclusions, it is necessary to approximate the solution of the background problem (<ref>) and
    the adjoint problem (<ref>). Moreover, when considering synthetic data u_meas, we must be able to compute the solution to the perturbed problem
    (<ref>) in presence of the exact inclusion. We rely on the Galerkin finite element method for the numerical approximation of these problems.
    The one-shot procedure makes the reconstruction algorithm very efficient, only requiring the solution of an adjoint problem for each acquired measurement over
    the time interval, without entailing any iterative (e.g. descent) method for numerical optimization.
    
    
    
     §.§ Finite Element approximation of  initial and boundary value problems
    
    
    The background problem (<ref>) can be cast in weak form as follows
    
    ∀ t∈(0,T), find u(t) ∈ V = H^1(Ω) such that u(0) = u_0 and
        ∫_Ω u_t(t)v + ∫_Ω k_1 ∇ u(t)·∇ v + ∫_Ωf(u(t))v = 0,     ∀  v ∈ V.
    
    By introducing a finite-dimensional subspace V_h of V, dim(V_h) = N_h < ∞, the Galerkin (semi-discretized in space) formulation of problem
    (<ref>) reads
    
    ∀  t∈(0,T), find u_h(t) ∈ V_h such that u_h(0) = u_h,0 and
        ((u_h)_t(t),v_h) + b(u_h(t),v_h) + F(u_h(t),v_h) = 0,     ∀  v_h ∈ V_h,
    
    where (·,·) is the inner product in L^2(Ω), b(u,v) = ∫_Ωk_1 ∇ u ·∇ v, F(u,v)=∫_Ω f(u)v, f is defined as in (<ref>) and u_h,0 is the H^1-projection of u_0 onto V_h.
    
    
    	
    
    
    
    
    
    	
    
    
    
    
    To obtain a full discretization of the problem, we introduce a finite difference approximation in time. According to the strategy reported in <cit.>, <cit.>,
    we rely on a semi-implicit scheme which allows an efficient treatment of the nonlinear terms.
    Let us consider an uniform partition {t^n}_n = 0^N of the time interval [0,T] of step τ=T/N s.t. t^0 = 0,  t^N=T. Then, the fully discrete
    formulation of (<ref>) is given by
    
    ∀ n = 0, … N-1, find u_h^n+1∈ V_h such that u_h^0 = u_0,h and
        (u^n+1_h,v_h) - (u^n_h,v_h) + τ b(u^n+1_h,v_h) + τ F(u^n_h,v_h) = 0,     ∀  v_h ∈ V_h.
     With the same discretization strategy one may describe a numerical scheme for the approximate solution of the perturbed problem, using the weak form reported in (<ref>) and introducing the forms
    
        b_ε(u,v) = ∫_Ω k_ε∇ u ·∇ v,     
        	F_ε(u,v) = ∫_Ωχ_Ω∖ω_ε f(u)v.
    
    The adjoint problem, instead, requires the introduction of the form dF(u,v;w) = ∫_Ω f'(w)uv, which is bilinear with respect to u and v. Thanks
    to the linearity of the adjoint problem, we can consider a fully implicit Crank-Nicolson scheme
    
    ∀  n = 0, … N-1, find w_h^n∈ V = H^1(Ω) such that w^N_h = 0 and
        (w^n+1_h,v_h) - (w^n_h,v_h) +τ/2( b(w^n+1_h,v_h) + b(w^n_h,v_h) + . 
           . dF(w^n+1_h,v_h;u^n+1_h) + dF(w^n_h,v_h;u^n_h) )= 
           τ/2( ∫_∂Ω (u_h^n+1 - u_meas(t^n+1))v_h + ∫_∂Ω (u_h^n - u_meas(t^n))v_h  ),   ∀  v_h ∈ V_h.
    
    The existence and the uniqueness of the solutions of the fully-discrete problems (<ref>) and (<ref>) follow by the well-posedness of the continuous problems, since V_h is a subspace of H^1(Ω).
    For further details on the stability and of the convergence of the proposed schemes we refer to <cit.>, <cit.> and <cit.>.
    
    
    The numerical setup for the simulation is represented in Figure <ref>. We consider an idealized geometry of the left ventricle (which has been object of several studies, see e.g. <cit.>,
    <cit.>), and define a tetrahedral tesselation 𝒯_h of the domain. The discrete space V_h is the P1-Finite Element space over 𝒯_h, i.e. the space of the continuous functions
    over Ω which are linear polynomials when restricted on each element T ∈𝒯_h. The mesh we use for all the reported results consists of 24924 tetrahedric elements and N_h=5639 nodes.
    We report also the anisotropic structure considered in all the recontruction tests, according to <cit.> and <cit.>. The conductivity matrix 𝕂_0 for the monodomain equation is given by 𝕂_0(x) = 𝕂^e(x)(𝕂^e(x)+ 𝕂^i(x))^-1𝕂^i(x), where both 𝕂^i and 𝕂^e are orthotropic tensors with three constant positive real eigenvalues, namely
    
        𝕂^e(x)    = k_f^e e⃗_⃗f⃗(x)⊗e⃗_⃗f⃗(x) + k_t^e e⃗_⃗t⃗(x)⊗e⃗_⃗t⃗(x) + k_r^e e⃗_⃗r⃗(x)⊗e⃗_⃗r⃗(x)
        
        𝕂^i(x)    = k_f^i e⃗_⃗f⃗(x)⊗e⃗_⃗f⃗(x) + k_t^i e⃗_⃗t⃗(x)⊗e⃗_⃗t⃗(x) + k_r^i e⃗_⃗r⃗(x)⊗e⃗_⃗r⃗(x)
    
    The eigenvectors e⃗_⃗f⃗, e⃗_⃗t⃗ and e⃗_⃗r⃗ are associated to the three principal directions of conductivity in the heart tissue: respectively, the fiber centerline, the tangent direction to the heart sheets and the transmural direction (normal to the sheets).
    
    
    
    For the direct problem simulations, we consider the formulation reported in (<ref>), specifying realistic values for the parameters C_m and ν. We have rescaled the values of the coefficients u_1, u_2, u_3 and A^2 in order to simulate the electric potential in the adimensional range [0,1]. The rescaling is performed by the transformation ũ = (α+u)/β, where α = 0.085 mV and β = 0.125 mV, whereas for the sake of simplicity we will still denote by u the rescaled variable ũ. We consider the initial datum u_0 to be positive on a band of the endocardium, representing the initial stimulus provided by the heart conducting system. The most important parameters, considered in accordance with <cit.>, <cit.>, are reported in Table <ref>.
    
    
    
    In Figure <ref> we report the solution of the discrete background problem (<ref>) at different time instants, comparing the isotropic and the anisotropic cases.
    
    
     §.§ Reconstruction of small inclusions
    
    We now tackle the problem of reconstructing the position of a small inhomogeneity using the knowledge of the electric potential of the tissue on a portion Γ
    of the boundary of the domain. In particular, we assume that u_meas is known on the endocardium, i.e. the inner surface of the heart cavity. We generate
    synthetic data on a more refined mesh and test the effectiveness of Algorithm 1 in the reconstruction of a small spherical inclusion in different positions.
    In Figure <ref> 
    we report the value assumed by the topological gradient, and superimpose the exact inclusions: we observe a negative region in proximity of the position of the real inclusion.
    The algorithm precisely identifies the region where the inclusion is present, whereas the minimum may in general be found along the endocardium also when the center of the real inclusion is not located on the heart surface.
    Nevertheless, due to the thinness of the domain the reconstructed position is close to the real one.
    
    This slight loss in accuracy seems to be an intrinsic limit of the topological gradient strategy applied to the considered problem. We point out that the reconstruction is performed by relying on a single measurement
    acquired on the boundary, which is indeed a constraint imposed by the physical problem. Hence, all the techniques relying on the introduction of several measurements to increase the quality of the reconstruction are impracticable.
    A different strategy, as proposed in several works addressing the steady-state case, may consist in introducing a modification to the cost functional J. In <cit.> and many related works the authors
    introduce a cost functional inherited from imaging techniques, whereas in <cit.>, <cit.> different strategies involving the Kohn-Vogelius functional or similar ones are explored. Nevertheless,
    the nonlinearity of the direct problem considered in this work prevents the possibility to apply these techniques, since the analitycal expressions of the fundamental solution, single and double layer potentials would not be available.
    
    
    
     §.§ Reconstruction in presence of experimental noise
    
    We test the stability of the algorithm in presence of experimental noise on the measured data u_meas. We consider different noise levels, according to the formula
    
        u_meas(x,t) = u_meas(x,t) + p η(x,t),
    
    where η(x,t), for each point x and instant t, is a Gaussian random variable with zero mean and standard deviation equal to u_3-u_1, whereas p ∈ [0,1] is the noise level.
    In Figure <ref> the results of the reconstruction with different noise levels are compared. The algorithm shows to be highly stable with respect to high rates of noise, with increasing accuracy as the noise level reduces.
    
    
     §.§ Reconstruction from partial discrete data
    
    A further test case we have performed deals with the reconstruction of the position of small inclusions starting from the knowledge of partial data. We are interested in assessing the effectiveness of our algorithm when the electric potential is measured only on a discrete set of points on the endocardium, possibly simulating the procedure of intracavitary electric measurements. Figure <ref> shows that the algorithm is able to detect the region where the small ischemia is located from the knowledge of the potential on N_p = 246, 61, 15 different points. The
    position of the reconstructed inclusion is slightly affected by the reduction of sampling points; nevertheless, reliable reconstructions can be obtained even with a very small number of points.
    
    
    	
    For the same purpose, we have tested the capability of the reconstruction procedure to avoid false positives: the algorithm is able to distinguish the presence of a real ischemia from the case where no ischemic region is present, also in the case where the data are recovered only at a finite set of points, and are affected by noise. We compare the value of the cost functional and of the minimum of the topological gradient obtained through Algorithm <ref> on data generated when (i) a small ischemia is present in the tissue or (ii) no inclusion is considered. The measurement is performed on a set of N_p = 100 points and is affected by different noise levels. The results are reported in Table <ref>.
    
    We observe that the presence of a small noise on the measured data causes a great increase of the cost functional J: with 5% noise, e.g., the value of J is two orders of magnitude greater than the value that J assumes in presence of a small inclusion without noise.
    Nevertheless, the topological gradient G allows to distinguish the false positive cases, since (at least in case of a small noise level) the value attained by its minimum in presence of a small inclusion is considerably lower than the random oscillations of G due to the noise.
    
    
     §.§ Reconstruction of larger inclusions
    
    We finally assess the performance of Algorithm 1, developed for the reconstruction of small inclusions well separated from the boundary, in detecting the position of extended inclusions, which may be of greater interest in view of the problem of detecting ischemic regions.
    Indeed, total occlusion of a major coronary artery generally causes the entire thickness of the ventricular wall to become ischemic (transmural ischemia) or, alternatively, a significant ischemia only in the endocardium, that is, the inner layer of the myocardium (subendocardial ischemia). See, e.g., <cit.> for a detailed investigation of the interaction between the presence of moderate or severe subendocardial ischemic regions and the anisotropic structure of the cardiac muscle.
    
    
    The most important assumption on which our one-shot procedure relies is that the variation of the cost functional from the value J(0) attained in the background case to the value related to the exact inclusion can be correctly described by the first order term of its asymptotic expansion, the topological gradient G. Removing the hypothesis of the small size, we cannot rigorously assess the accuracy of the algorithm, however  it still allows us to identify the location of the ischemic region.
    
    The results reported in Figure <ref> show that in presence of a  inclusion of larger size (and not even separated from the boundary), the minimum of the topological gradient is found to be close to the position of the inclusion, and attains lower values with respect to the previously reported cases.
    
    Moreover, in Figure <ref> we also assess the stability of the reconstruction with respect to the presence of noisy data and partial measurements, as done in the case of small inclusions.
    
    
    § CONCLUSIONS AND PERSPECTIVES
    
    A rigorous theoretical analysis of the inverse problem of detecting inhomogeneities in the monodomain equations has allowed us to set up a numerical reconstruction procedure, aiming at the detection of ischemic
    regions in the myocardic tissue from a single measurement of the endocardial potential. The identification is made possible by evaluating the topological gradient  of a quadratic cost functional, requiring
    the solution of two initial and boundary value problems, the background problem and the adjoint one. Numerical results are encouraging and allow to estimate the position of the inclusion,
    although the identified inhomogeneity is nearly always detected on the boundary where the measurement is acquired. Nevertheless, provided a single measurement can be used for the sake of  identification,
    and a  one-shot procedure is performed, the obtained results give useful insights.
    
    Many issues are still open. Concerning the mathematical model, an even more interesting case would be the one involving the heart-torso coupling is considered, so that more realistic (and noninvasive)
    body surface measurements can be employed. Setting and analyzing the inverse problem in this context represents the natural continuation  of the present work.
    To close the gap between the rigorous mathematical setting and the practice, the two assumptions made in this work about the size of the inclusion and its distance from the boundary should be relaxed.
    Numerical results shown in Section <ref>  provide a first insight on the detection of inclusions with larger size, as those corresponding to transmural or subendocardial ischemias.
    From a mathematical standpoint, this problem is still open. Also in the case of a linear direct problem, very few results can be found in literature, see, e.g. <cit.>.
    Estimating the size of the inclusion is another open question in the case of parabolic PDEs, also for linear equations. The case of multiple inclusions, addressed in <cit.>
    for a stationary nonlinear problem, could also be considered. Last, but not least, the topological optimization framework addressed in this paper could also be combined with an iterative algorithm,
    such as the level set method, or with the solution of a successive shape optimization problem, to achieve a full reconstruction both of the dimension and the shape of the inclusion.
    
    
    
    
    
    
    § APPENDIX - PROOF OF PROPOSITION <REF>
    
    Setting Z(t) = Φ(T-t),   t ∈ (0,T), we get an equivalent problem to (<ref>)
        Z_t - k_0Δ Z + f^'(u)Z  = 0,    in  Ω× (0,T),
           ∂ Z/∂ n  = u^ε - u,     on   ∂Ω× (0,T), 
            Z(0) = 0,     in  Ω.
        			2truecm
    
    We will prove that Z ∈ L^2(0,T;H^3(K)) ↪ L^1(0,T;W^1,∞(K)). To this aim we need to derive some energy estimates.
    Multiplying the first equation in (<ref>) by Z, an application of Young's inequality leads to
    
        1/2d/dtZ(t)^2_L^2(Ω)  + k_0/2∇ Z_L^2(Ω)^2
        ≤ C (Z_L^2(Ω)^2 + u^ε -u^2_L^2(∂Ω)),
    
    where C= C(k_0, M_2, Ω) >0. An application of Gronwall's lemma gives
    
        Z(t)^2_L^2(Ω))≤ C u^ε -u^2_L^2(0,t; L^2(∂Ω)),   ∀  t ∈ [0,T],
    
    so that
    
        Z^2_L^∞(0,t; L^2(Ω))≤ Cu^ε -u^2_L^2(0,t; L^2(∂Ω)),   ∀  t ∈ [0,T].
    
    Instead, integrating (<ref>) in time over [0,t] we get
    
        ∫_0^t∇ Z_L^2(Ω)^2
        ≤ C ( ∫_0^tZ_L^2(Ω)^2 + ∫_0^tu^ε -u^2_L^2(∂Ω) )
    
    and finally
    
        ∇ Z_L^2(0,t;L^2(Ω))^2
        ≤ Cu^ε -u^2_L^2(0,t;L^2(∂Ω)),   ∀  t ∈ [0,T],
    
    where C is a positive constant depending on k_0, M_2, Ω, T.
    We remark that, by standard regularity results, Z is smooth on E × [0,T], for any compact E⊂Ω.
    
    Consider now two compact sets K_1 and K_2 such that
    
        K⊂ K_2 ⊂ K_1 ⊂Ω,    d(k_0, ∂Ω) ≥ d_1 >0.
    
    It is possible to construct two functions ξ_1,  ξ_2 and two constants b_1,   b_2 satisfying
    
        ξ_i ∈ C^2(Ω),    0≤ξ_i ≤ 1,   ξ_i(x) = 1   ∀  x ∈ K_i,   ξ_i(x) = 0
          ∀ x ∈ B_i   i=1,2,
    
        B_i = { x ∈Ω  :   d(x, ∂Ω) ≤ b_i},    0 < b_1 < b_2 < d_1,   
        K ⊂⊂ Supp ξ_2 ⊂⊂ K_1 ⊂ Supp ξ_1 ⊂Ω.
    
    Let us multiply the first equation of (<ref>) by -Δ Z. Then the following holds
    
        d/dt (1/2|∇ Z|^2) + k_0(Δ Z)^2 - f^'(u)ZΔ Z =  div (Z_t ∇ Z).
    
    Multiplying (<ref>) by ξ_1, integrating on Ω× (0,T) and using the definitions of Z we get
    
        ∫_Ω (1/2|∇ Z(T)|^2 )ξ_1 + k_0∫_0^T∫_Ω(Δ Z)^2ξ_1
        = ∫_0^T∫_Ω f^'(u)ZΔ Z ξ_1  - ∫_0^T∫_Ω Z_t ∇ Z ·∇ξ_1.
    
    Combining (<ref>) and the first equation in (<ref>), applying Young's inequality and taking into account (<ref>) and the fact
    that 0 ≤ξ≤ 1, we obtain
    
        ∫_Ω (|∇ Z(T)|^2 )ξ_1 + k_0∫_0^T∫_Ω(Δ Z)^2ξ_1
        ≤ 2M_2∫_0^T∫_Ω Z^2- 2∫_0^T∫_Ω (k_0Δ Z - f^'(u)Z) ∇ Z ·∇ξ_1.
    
    Integrating by parts the term ∫_0^T∫_ΩΔ Z ∇ Z ·∇ξ_1, we can easily deduce
    
    
    
    
    
    
    
        ∫_Ω (|∇  Z(T)|^2 )ξ_1 + ∫_0^T∫_Ω(Δ  Z)^2ξ_1
        ≤ C ( Z_L^2(0,T;L^2(Ω))^2 + ∇ Z^2_L^2(0,T;L^2(Ω))),
    
    where C is a positive constant depending on M_2, k_0, ξ_1.
    Hence, since ξ_1 = 1 in k_0, we get
    
        Δ  Z ^2_L^2(0,T;L^2(K_1))≤ C (  Z_L^2(0,T;L^2(Ω))^2 + ∇ Z^2_L^2(0,T;L^2(Ω))).
    
    Observe that, replacing T by t ∈ (0,T] in (<ref>), we deduce also
    
        ∇  Z_L^∞(0,T;L^2(K_1))≤ C ( Z_L^2(0,T;L^2(Ω))^2 + ∇ Z^2_L^2(0,T;L^2(Ω))).
    
    Combining (<ref>), (<ref>) and (<ref>), we obtain
    
        Z ^2_L^2(0,T;H^2(K_1))≤ C u^ε - u_L^2(0,T;L^2(∂Ω))^2,
    
    where C is a positive constant depending on k_0, M_2, Ω, T, ξ_1.
    
    On account of the first equation in (<ref>) and the previous estimates, we get
    
        Z_t_L^2(0,T;L^2(K_1))^2  ≤ C ( Z_L^2(0,T;L^2(Ω))^2
        + ∇  Z_L^2(0,T;L^2(Ω))^2) ≤  C u^ε - u_L^2(0,T;L^2(∂Ω))^2,
    
    where C is a positive constant depending on k_0, M_2, Ω, T, ξ_1.
    
    Now, let us multiply the first equation of (<ref>) by  -Δ  Z_t. We obtain
    
        -  Z_t Δ  Z_t + k_0/2d/dt(Δ Z)^2 - f^'(u) ZΔ  Z_t = 0.
    
    Multiplying the previous equation by ξ_2 and integrating on Ω× (0,T), then a
    suitable integration by parts in space implies
    
        ∫_0^T∫_Ω|∇ Z_t|^2ξ_2 + k_0/2∫_0^T∫_Ωd/dt(Δ Z)^2ξ_2
        + ∫_0^T∫_Ωξ_2 Z f^''(u)∇ u ·∇ Z_t + ∫_0^T∫_Ωξ_2f^'(u)∇  Z·∇ Z_t
    
        = ∫_0^T∫_Ω div  ( 1/2∇(( Z_t)^2) ) ξ_2
        + ∫_0^T∫_Ω div  (∇(f^'(u) Z  Z_t) -  Z_t∇(f^'(u) Z))ξ_2.
    
    Integrating by parts the second term of the left-hand side and by parts in space the terms in the right-hand side, by an application of Young's
    inequality we finally get
    
        ∫_0^T∫_K_2|∇ Z_t|^2 ≤∫_0^T∫_Ω|∇ Z_t|^2ξ_2 ≤
        C ( ∫_0^T∫_Ω| Z|^2  + ∫_0^T∫_Ω  |∇ Z|^2 + ∫_0^T∫_K_1 ( Z_t)^2 ),
    
    where the constant C>0 depends on ξ_2, M_2.
    A combination with (<ref>), (<ref>), (<ref>) gives
    
        ∇ Z_t_L^2(0,T;L^2(K_2))^2 ≤
        C u^ε - u _L^2(0,T;L^2(∂Ω))^2,
    
    where the constant C>0 depends on k_0, M_2, Ω, T, ξ_1, ξ_2.
    In order to prove the desired regularity, we need to take into account also the third-order derivatives, in particular the operator ∇Δ  Z.
    Observe that from the first equation in (<ref>) we get
    
        ∇Δ  Z = 1/k_0 ( ∇ Z_t +  Z f^''(u)∇ u + f^'(u)∇  Z  ).
    
    Hence, we can conclude
    
        ∇Δ Z_L^2(0,T;L^2(K_2))^2 ≤
        C u^ε - u _L^2(0,T;L^2(∂Ω))^2,
    
    where C is a positive constant depending on k_0, 1/k_0, M_2, Ω, T, ξ_1, ξ_2.
    
    Recalling (<ref>) and the fact that K ⊂ K_2 ⊂ K_1, standard regularity results imply
    
        Z_L^2(0,T;H^3(K))^2 ≤
        C u^ε - u _L^2(0,T;L^2(∂Ω))^2.
    
    Finally, from (<ref>) and (<ref>), by Sobolev immersion theorems, we get
    
        Z_L^1(0,T;W^1,∞(K))^2  ≤ C(T)  Z_L^2(0,T;W^1,∞(K))^2 ≤
        C u^ε - u _L^2(0,T;L^2(∂Ω))^2,
    
    where C is a positive constant depending on k_0, 1/k_0, M_2, Ω, T, ξ_1, ξ_2.
    
    Recalling the relation between Z and Φ we get (<ref>).
    
    
    § ACKNOWLEDGMENTS
    
    E. Beretta, C. Cavaterra, M.C. Cerutti and L. Ratti thank the New York University in Abu Dhabi
    for its kind hospitality that permitted a further development of the present research.
    The work of C. Cavaterra was supported by the FP7-IDEAS-ERC-StG 256872 (EntroPhase)
    and by GNAMPA (Gruppo Nazionale per l'Analisi Matematica, la Probabilità
    e le loro Applicazioni).
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    99Alvarez2012 D. Álvarez, F. Alonso-Atienza, J.L. Rojo-Álvarez,  A. Garcia-Alberola, and M. Moscoso,
    Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: a model study,  Math. Computer Modeling
    55 (2012) 1770–1781.
    
    Ammari2013
    H. Ammari, P. Garapon, F.  Jouve, H. Kang, M. Lim, and S. Yu,
     A new optimal control approach for the reconstruction of extended inclusions, SIAM J. Control Optim.  51(2) (2013)  1372–1394.
    
    ammari2012stability H. Ammari, J. Garnier, V. Jugnon, and H. Kang, Stability and resolution analysis for a topological derivative based imaging functional, SIAM J. Control Optim.
    50 (2012) 48–76.
    
    Ammari2005 H. Ammari, E. Iakovleva, H. Kang, and K. Kim, Direct algorithms for thermal imaging of small inclusions, Multiscale Model. Simul.  4 (4) (2005) 1116–1136.
    
    AK H. Ammari and H. Kang,  Reconstruction of small inhomogeneities from boundary measurements, Lectures Notes in Mathematics Series 1846 Springer, 2004.
    
    Bendahmane2006 M. Bendahmane and   K.H. Karlsen, Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue,
    Netw. Heterog. Media 1 (2006) 185–218.
    
    BCMP E. Beretta, M. C. Cerutti, A. Manzoni, and D. Pierotti, An asymptotic formula for boundary potential perturbations
    in a semilinear elliptic equation related to cardiac electrophysiology, Math. Mod. Meth. Applied Sciences 26 (2016) 645–670.
    
    BerettaManzoniRatti E. Beretta, A. Manzoni, and L. Ratti, A reconstruction algorithm based on topological gradient for an inverse problem related to a semilinear elliptic boundary value problem, Inverse Problems, accepted for publication  (2017).
    
    Gerbeau2008 M. Boulakia,  M.A Fernández, J.F. Gerbeau, and N. Zemzemi,  A coupled system of PDEs and ODEs arising in electrocardiograms modeling,
    Applied Math. Res. Exp. 2008 (2008), doi 10.1093/amrx/abn002.
    
    Bourgault2009 Y. Bourgault, Y.  Coudière, and C. Pierre, Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology,
    Nonlinear Anal. Real World Appl. 10 (2009) 458–482.
    
    Burger2010 M. Burger, K. A. Mardal, and B. F. Nielsen, Stability analysis of the inverse transmembrane potential problem in electrocardiography, Inverse Problems  26 (2010) 105012.
    
    CV Y. Capdeboscq and M. Vogelius, A general representation formula for boundary voltage perturbations caused by internal conductivity
    inhomogeneities of low order fraction, Math. Modelling and Num. Analysis 37 (2003) 159–173.
    
    CMV D.J. Cedio-Fengya, S. Moskow, and M. Vogelius, Identification of conductivity imperfections of small diameter by boundary measurements.
    Continuous dependence and computational reconstruction. Inverse Problems 14 (1998), no. 3, 553595.
    
    art:cmm S. Chaabane, M. Masmoudi, and H. Meftahi, Topological and shape gradient strategy for solving geometrical inverse problems,
    J. Math. Anal. Appl. 400 (2013) 724–742.
    
     Chavez2015
     C.E. Chávez, N. Zemzemi, Y. Coudière, F.  Alonso-Atienza, and D.  Álvarez, Inverse Problem of Electrocardiography: Estimating the Location of Cardiac Ischemia in a 3D Realistic Geometry, in Functional Imaging and Modeling of the Heart: 8th International Conference, FIMH 2015, Maastricht, The Netherlands, June 25-27, 2015. Proceedings, H. van Assen, P. Bovendeerd and T. Delhaas (Eds.), Springer International Publishing (2015)
    
    
    colli2004parallel P. Colli Franzone and L. Pavarino, A parallel solver for reaction–diffusion systems in computational electrocardiology, Math. Mod. Meth. Applied Sciences
    14(6) (2004) 883–911.
    
    Pavarino2010 P. Colli Franzone, L. Pavarino, and S. Scacchi, Dynamical effects of myocardial ischemia in anisotropic cardiac models in three dimensions, Math. Mod. Meth. Applied Sciences 17(12) (2007) 1965–2008.
    
    
    Pavarino2014book P. Colli Franzone, L. Pavarino, and S. Scacchi, Mathematical Cardiac Electrophysiology, Modeling, Simulation and Applications (MS&A) Series,
    13  Springer-Verlag Italia, Milano, 2014.
    
    ColliFranzone1978 P. Colli Franzone, B. Taccardi, and C. Viganotti, An approach to inverse calculation of epicardial potentials from body surface maps,
    Adv. Cardiol. 21 (1978) 50–54.
    
    DiCristoVessella2010 M. Di Cristo and S. Vessella, Stable determination of the discontinuous conductivity coefficient of a parabolic equation, SIAM J. Math. Anal.  42(1) (2010) 183–217.
    
    D L. Dung, Remarks on Hölder continuity for parabolic equations and convergence to global attractors, Nonlinear Analysis 41 (2000) 921–941.
    
    Isakov1 A. Elayyan and V. Isakov, On uniqueness of recovery of the discontinuous conductivity coefficient of a parabolic equation, SIAM J. Math. Anal.  28(2) (1997) 49–59
    
    fernandez2010decoupled M. Fernández and N. Zemzemi, Decoupled time-marching schemes in computational cardiac electrophysiology and ECG numerical simulation, Math. Biosci.
    226(1) (2010) 58–75.
    
    FV A. Friedman and M. Vogelius, Identification of small inhomogeneities of extreme conductivity by boundary measurements: a theorem on continuous dependence, Arch. Rational Mech. Anal. 105 (1984), 299-326.
    
    gerbeau2015reduced J.F. Gerbeau, D. Lombardi,  and E. Schenone, Reduced order model in cardiac electrophysiology with approximated Lax pairs, Adv. Comput. Math.
    41(5) (2015) 1103–1130.
    
    GT D. Gilbarg and N.S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer Verlag, Berlin, 1983.
    
    Isakov2 V. Isakov, K. Kim, and G. Nakamura,   Reconstruction of an unknown inclusion by thermography, Ann. Scuola Norm. Sup. Pisa Cl. Sci.  (5) Vol. IX (2010) 725–758.
    
    LSU O.A. Ladyzenskaja, V.A. Solonnikov, and N.N. Ural'ceva, Linear and Quasi-linear Equations of Parabolic Type, AMS Transl. Monographs 23
    AMS, Providence, 1968.
    
    L A. Lunardi, Analytic semigroups and optimal regularity in parabolic problems,  Birkhäuser, Basel, 1995.
    
    LN M. Lysaker and B.F. Nielsen,  Towards a level set framework for infarction modeling: an inverse problem, Int. J. Numer. Anal. Model.
    3 (2006) 377–394.
    
    NielsenEtAl20072 B.F. Nielsen, X. Cai, and M. Lykaser, On the possibility for computing the transmembrane potential in the heart with a one shot method: an inverse problem,
    Math. Biosciences 210 (2007) 523–553.
    
    Nielsen2007  B.F. Nielsen, M. Lykaser, and A. Tveito,  On the use of the resting potential and level set methods for identifying ischemic heart disease: An inverse problem,
    J. Comput. Phys. 220 (2007) 772–790.
    
    N L. Nirenberg, On elliptic partial differential equations, Ann. Scuola Norm. Sup. Pisa 13 (1959) 115–162.
    
    P C.V. Pao, Nonlinear parabolic and elliptic equations,  Plenum Press, New York, 1992.
    
    park2012topological  W.K. Park, Topological derivative strategy for one-step iteration imaging of arbitrary shaped thin, curve-like electromagnetic inclusions, J Comput. Phys.
    231(4) (2012) 1426–1439.
    
    MacLeod2010 A.J. Pullan,  L.K. Cheng, M.P. Nash, A. Ghodrati, R. MacLeod, and D.H. Brooks, The inverse problem of electrocardiography, Comprehensive Electrocardiology,
    299–344,  P.W. Macfarlane and A. van Oosterom and O. Pahlm and P. Kligfield and  M. Janse and J. Camm editors, Springer London, 2010.
    
    Quarteroni2016 A. Quarteroni, T. Lassila, S. Rossi, and R. Ruiz-Baier, Integrated Heart – Coupling multiscale and multiphysics models for the simulation of the cardiac function, Comput. Methods Appl. Mech. Engrg. 314 (2017) 345–407.
    
    R J.C. Robinson, Infinite-Dimensional Dynamical Systems, Cambridge texts in applied mathematics, Cambridge University Press, Cambridge, 2001.
    
    rossi2014thermodynamically S. Rossi, T. Lassila, R. Ruiz-Baier, A. Sequeira,  and A. Quarteroni, Thermodynamically consistent orthotropic activation model capturing ventricular systolic wall thickening in cardiac electromechanics, Eur. J. Mech. A-Solid.
    48 (2014) 129–142.
    
    sanfelici2002convergence S. Sanfelici, Convergence of the Galerkin approximation of a degenerate evolution problem in electrocardiology, Numer. Methods Partial Differential Equations
    18(2) (2002) 218–240.
    
    SLCNMT J. Sundnes, G.T. Lines, X. Cai, B.F. Nielsen, K.A. Mardal, and A. Tveito, Computing the electrical activity in the heart,
    Monographs in Computational Science and Engineering Series, 1, Springer, (2006).
