Well-posedness for a modified bidomain model describing bioelectric activity in damaged heart tissues

We prove the existence and the uniqueness of a solution for a modified bidomain model, describing the electrical behaviour of the cardiac tissue in pathological situations. The leading idea is to reduce the problem to an abstract parabolic setting, which requires to introduce several auxiliary differential systems and a non-standard bilinear form. The main difficulties are due to the degeneracy of the bidomain system and to its non-standard coupling with a diffusion equation, accounting for the presence of the pathological zone in the heart tissue.


Introduction
In this paper, we are interested in studying a modified version of the famous bidomain model (see, e.g., [10,19,20] and the references therein; see, also, the references quoted in [11,Introduction]), which is one the most well-known mathematical models in cardiac electrophysiology. This is a topic of major interest in biomedical research.
In the classical bidomain model, at a macroscopic scale, the electric activity of the heart is governed by a system of two degenerate reaction-diffusion partial differential equations for the averaged intra-cellular and, respectively, extra-cellular electric potentials, along with the transmembrane potential, coupled in a nonlinear manner to ordinary differential equations describing the dynamics of the ion channels. The well-posedness of the bidomain model has been studied, for different nonlinear 1 ionic models and by using different techniques, by several authors (see, for instance, [3,8,11,14,15,18,21,23,24]). The bidomain model is suitable for describing the propagation of the action potential in a perfectly healthy cardiac tissue, but it is no longer valid (even if one tries to ad-hoc modify some of its relevant modeling parameters) in pathological situations. Models, taking into account the presence in the cardiac tissue of damaged zones, called diffusive inclusions and assumed to be passive electrical conductors, were proposed in [6,12,13,14]. From a mathematical point of view, such models consist in a bidomain system coupled with a diffusion equation. More precisely, one has a degenerate reaction-diffusion system of partial differential equations modeling the intra-cellular and, respectively, the extra-cellular electric potentials of the healthy cardiac tissue, coupled with an elliptic equation for the passive regions and with an ordinary differential equation describing the cellular membrane dynamics. We point out that in all the above mentioned papers a perfect electrical coupling between the healthy part of the heart and the damaged tissue was assumed. More general conditions for the heart-torso coupling were proposed in [7] and investigated through numerical simulations in [4,5,26], in order to take into account the possible capacitive and resistive effects of the pericardium. However, up to our knowledge, there are no rigorous proofs in the literature covering this setting. We investigate these more general conditions in the context of the bidomain model with diffusive inclusions, where the appropriate interface behaviour, up to our knowledge, is still not well understood. The goal of the present paper is to study the well-posedness of such a modified bidomain model. We include the structural defects of the heart tissue in this model by coupling a standard bidomain system in the healthy zone with a diffusion equation posed in the damaged part of the heart, through non-standard conditions (see equations (2.15)-(2.20)). More precisely, for the intra-cellular potential we assume no flux condition on the interface between the two zones (see (2.18)), while the extra-cellular potential is coupled with the electrical potential of the damaged zone through imperfect transmission conditions, involving the resistive and the capacitive properties of the interface (see (2.19) and (2.20)). In order to describe the dynamic of the membrane, one can use a physiological ionic model or a phenomenological one (see, for instance, [11]). In this paper, the dynamic of the gating variable modeling the ionic transport through the cell membrane is described with the aid of a Hodgkin-Huxley type formalism (see (2.10)-(2.13)). Our analysis covers also the modified Mitchell-Schaeffer formalism proposed in [13] (see Remark 2.3). We point out again that our mathematical model generalizes the modified bidomain model with diffusive inclusions and perfect transmission conditions considered in [13,14], the original model being recovered by suitably rearranging the parameters appearing in equation (2.20). We believe that further numerical simulations have to be carried out in order to validate the relevance of such transmission conditions also from the point of view of possible biological applications.
The mathematical problem we address here is rather non-standard and, up to our knowledge, the proof of its well-posedness is new in the literature and generates difficulties due to the degeneracy of the bidomain system and to its special coupling with the diffusion equation. Our main result is contained in Theorems 3.4 and 3.6, where the leading idea is to reduce the problem to an abstract parabolic setting (see [9,22]). This requires to introduce several auxiliary differential systems and a non-standard bilinear form (see Proposition 3.3). The problem proposed here can be seen as a mesoscopic model which will be analyzed in the homogenization limit in a forthcoming article (see [2]).
The paper is organized as follows: in Section 2, we introduce the mathematical description of our modified bidomain model, together with its geometrical and functional setting. In Section 3, we state and prove our main result.

The model
Let Ω be an open connected bounded subset of R N ; we assume that ∂Ω is of class C ∞ , though this assumption can be weakened. Moreover, for T > 0, we set Ω T = Ω × (0, T ). We assume that where Ω D and Ω B are two disjoint open subsets of Ω and Γ = ∂Ω D ∩ Ω = ∂Ω B ∩ Ω. The domain Ω is occupied by the cardiac tissue, Ω B represents the healthy part of the heart tissue, modeled with the aid of a standard bidomain system, Ω D represents the diffusive region, accounting for the damaged part of the heart, and Γ is the common boundary of these two regions, assumed to be Lipschitz. From a geometrical point of view, we assume that Ω B is connected, while Ω D might be connected or disconnected. Indeed, we will consider two different cases: in the first one (to which we will refer as the connected/disconnected case, see Fig.1 on the left), we will assume Ω D ⊂⊂ Ω and Ω D is made by a finite number of connected components. In this case, Γ = ∂Ω D and ∂Ω B ∩ ∂Ω = ∅. In the second case (to which we will refer as the connected/connected case, see Fig.1 on the right), we will assume that both Ω D and Ω B are connected, with ∂Ω B ∩∂Ω = ∅ and ∂Ω D ∩ ∂Ω = ∅. Finally, let ν denote the normal unit vector to Γ pointing into Ω B . In the following, by γ we shall denote a strictly positive constant, which may depend on the geometry and on the other parameters of the problem; γ may vary from line to line.
where H is endowed with the scalar product (w, r), (w, s) H := here, α > 0 will be the constant appearing later in (2.20) and W is endowed with the scalar product (w, r), (w, s) W := We denote by (·, ·) 1/2 the standard scalar product on H 1/2 (Γ ). Moreover, we define the space endowed with the norm We recall that ∂Ω B ∩ ∂Ω is always non-empty, while ∂Ω D can intersect or not the boundary of Ω, depending on the geometry. For W ∈ X 1 0 (Ω), we have the following Poincaré inequality (see [1,Proposition 2]): where [W] = W | Ω B −W | Ω D and the last term is not necessary in the connected/connected case. Therefore, an equivalent norm on X 1 0 (Ω) is given by W 2 again, the last term can be dropped in the connected/connected case.

2.3.
Position of the problem. Let α, β be strictly positive constants and in Ω, for suitable strictly positive constants γ 0 , γ 0 . The assumption that σ B 1 , σ B 2 , σ D are scalar functions is used only in Section 3. Removing this assumption is not trivial. If we want to consider general bounded and symmetric matrices satisfying  [3,4,14]).
Let us consider a locally Lipschitz continuous function g : R 2 → R, such that g(p, 1) ≥ 0 and g(p, 0) ≤ 0. The example we have in mind here is a function of the form where a, b : R → R are positive, bounded and Lipschitz functions. Notice that the form of g in (2.10) is classical in this framework (see, for instance, [24]) and that g is Lipschitz continuous with respect to p and affine with respect to q. Let I ion : R 2 → R be given by in Ω B , and p ∈ L 2 (Ω B T ). Consider the gating equation Notice that, by classical results, the previous problem admits a unique solution w p ∈ H 1 (0, T ; L ∞ (Ω B )) and, from our assumptions, in Ω B . This is a standard result for ODEs, taking into account that the spatial variable plays here only the role of a parameter (for similar results, see, for instance, [11,13,16]). Moreover, from the previous assumptions, we can prove that there exists a strictly positive constant γ such that due to the Lipschitz dependence of w p on p and to the bound 0 ≤ w p (x, t) ≤ 1 a.e.
in Ω B T . We give here a complete formulation of the problem we shall address in this paper. The operators div and ∇ act only with respect to the space variable x.
where w is the solution of the gating equation (2.12), (2.13), with p = u B 1 − u B 2 . Remark 2.1 (Biological interpretation). In the previous system of equations, the coefficients σ B 1 , σ B 2 and σ D are the conductivities of the two healthy phases and of the damaged one, respectively, while α and β are given parameters related to the capacitive and the resistive behaviour of the interface Γ . The functions f 1 and f 2 , appearing in (2.15) and (2.16), respectively, represent the internal and the external current stimulus. The solutions u B 1 and u B 2 are the intra and the extra-cellular potentials of the healthy zone, while u D is the electrical potential of the damaged zone. The function u B 2 − u B 1 is the so-called transmembrane potential. Finally, the variable w, called the gating variable, describes the ionic transport through the cell membrane. The terms g and I ion are nonlinear functions, modeling the membrane ionic currents. For simplicity, we consider only one gating variable, but our results hold true also for the case in which the gating variable is vector valued.  We consider here a Hodgkin-Huxley type model (see (2.10)-(2.11)), as in [2,11,24]. However, we point out that the results obtained in this paper are also valid for a regularized version of the Mitchell-Schaeffer model proposed in [13] (see, also, [12,14,18] (2.14) and moving the integral containing I ion to the right-hand side, we get T ) ) , (2.24) where γ and δ are positive constants, δ can be chosen smaller than min(σ B 1 , σ B 2 ), and we have also applied Poincaré inequality to u B 1 and u B 2 . By absorbing into the lefthand side the first two terms in the last line of (2.24) and using Gronwall inequality, from the previous estimate, we obtain Reasoning in a similar way as done for (2.24), i.e. by multiplying the first equation by υ B 1 , the second one by υ B 2 , the third one by υ D , subtracting the second equation from the first one, adding the third one, integrating by parts, using the remaining equation of the previous system moving the integral containing I ion to the right-hand side and using Hölder inequality, we get where, in the last inequality, we used (2.14). We can conclude by using Gronwall inequality.
Notice that, by setting V = u B 1 − u B 2 , U = u B 2 a.e. in Ω B , U = u D a.e. in Ω D , and denoting by [·] the jump across Γ of the quantity in the square brackets, i.e.

Well-posedness
We will consider our problem in an abstract setting and to this purpose we need first "to move" the source f 1 − f 2 from (2.27) to (2.26), (2.31), and (2.34). Then, we introduce a bilinear form a on W ×W such that the problem is reduced to the abstract scheme (3.47). This form is constructed with auxiliary functions which are obtained in Proposition 3.1 and in Remark 3.2. In Proposition 3.3, we prove the necessary properties of the form a. In Theorem 3.4 and Proposition 3.5, we prove existence of solutions to the problem with I ion ≡ 0. Finally, the full result is obtained in Theorem 3.6, where the complete problem will be treated as a nonlinear perturbation of this case (see, for instance, [9,13,16,22]). We start by considering, for a.e. t ∈ (0, T ), the following auxiliary problem: Clearly, problem (3.1)-(3.3) is classical and admits a unique solution u ∈ H 1 (0, T ; H 1 null (Ω B )). Moreover, we extend u inside Ω D by zero, so that it has a nonzero jump on Γ ; i.e., u ∈ X 1 0 (Ω). On Γ , let us define and consider the problem for (v, u) ∈ L 2 (0, T ; H 1 null (Ω B )) × L 2 (0, T ; X 1 0 (Ω)) given by The weak formulation of the previous problem is given by for every ϕ B ∈ L 2 (0, T ; H 1 null (Ω B )) ∩H 1 0 (0, T ; L 2 (Ω B )), ϕ D ∈ L 2 (0, T ; X 1 0 (Ω)), where, as before, [ϕ D ] = ϕ 1 D − ϕ 2 D and [ϕ D ] ∈ H 1 0 (0, T ; L 2 (Γ )). The weak formulation (3.14) shall be complemented with the initial conditions. Indeed, as it will be proved in Theorem 3.4, we have v ∈ C 0 ([0, T ]; L 2 (Ω B )) and [u] ∈ C 0 ([0, T ]; L 2 (Γ )). The next step is to define a suitable bilinear form on W , which is continuous and coercive. To this purpose, we need the following result.
Proposition 3.1. Let (w, r) ∈ W be assigned. Then, there exists a unique solution (3.17) [W] = r, on Γ ; Moreover, there exists a constant γ > 0, depending on σ B 1 , σ B 2 , σ D , and the geometry, such that Proof. Uniqueness for problem (3.15)-(3.19) is a straightforward consequence of its linearity. In order to prove that a solution does exist, we first consider the following auxiliary problem: in Ω B ; (3.21) (3.23) Clearly, the previous problem admits a unique solution W 1 ∈ X 1 0 (Ω). Moreover, there exists a constant γ > 0, depending on σ B 1 , σ B 2 and the geometry, such that Indeed, let us denote by r ∈ H 1 null (Ω B ) an extension of r from Γ to the whole Ω B , such that r H 1 (Ω B ) ≤ γ r H 1/2 (Γ ) , and set W r Therefore, by the standard energy inequality, we get which implies (3.25). Now, let us consider the second auxiliary problem for W 2 ∈ H 1 0 (Ω) given by Existence and uniqueness for the previous problem is guaranteed by [1,Lemma 5]; moreover, the weak formulation of (3.29)-(3.33) is given by which, replaced in (3.34), provides By taking ϕ = W 2 in (3.35), we get which, taking into account (3.25), implies with γ depending only on σ B 1 , σ B 2 , σ D and the geometry. Finally, setting W = W 1 + W 2 , it is easy to see that W ∈ X 1 0 (Ω) and satisfies (3.15)- (3.19) and also (3.20).
Remark 3.2. We point out that, taking r ∈ H 1 null (Ω B ) as in the proof of Proposition 3.1 and defining we get that W r ∈ H 1 0 (Ω) and satisfies The weak formulation of problem (3.38)-(3.41) is given by for every φ ∈ H 1 0 (Ω), which implies where W and W are the solutions of (3.15)- (3.19) corresponding to (w, r) and (w, s), respectively. Proof. Notice that the bilinear form a can be rewritten as which immediately proves that it is symmetric. Moreover, from (3.20), it easily follows that a is continuous. In order to prove that it is also coercive, we note that where, in the last inequality, we take into account that r = [W] (see (3.18)) and we use the Poincaré inequality (2.7) and the classical trace inequality, which assure that Proof. Let us denote by ·, · W * W the duality pairing between W and its dual space W * and define B ∈ W * as 15 for every (ϕ B , [ϕ D ]) ∈ W . By using the bilinear form a introduced in (3.45), we can consider the following abstract problem: in the distributional sense. By [25,Theorem 23.A], the problem (3.47) is well-posed and it is not difficult to see that its weak formulation is given by for every (ϕ B , [ϕ D ]) ∈ W and every ψ ∈ C ∞ 0 (0, T ), where ϕ D inside Ω B and Ω D is defined as the solution of (3.15)- (3.19), starting from ϕ B ∈ H 1 null (Ω B ) and [ϕ D ] ∈ H 1/2 0 (Γ, Ω). Clearly, (3.48) shall be complemented with the initial conditions. Notice that (3.48) formally coincides with (3.14); however, in (3.14) the test function ϕ D is a generic function belonging to L 2 (0, T ; X 1 0 (Ω)), while in the present case it is the solution of an assigned differential problem. Hence, in order to state that, actually, the two weak formulations are equivalent, we have to prove that we can replace the prescribed ϕ D in (3.48) with a generic test function belonging to X 1 0 (Ω). To this purpose, let us fix [ϕ D ] ∈ H 1/2 0 (Γ, Ω) and choose two generic functions ϕ 1 ∈ H 1 null (Ω B ) and (3.49) Hence, replacing (3.49) in (3.48), it follows that it is possible to take ϕ D ∈ X 1 0 (Ω) arbitrarily in (3.48) and, thus, such a weak formulation coincides with (3.14), once we take into account the density of product functions in L 2 (0, T ; H 1 null (Ω B )) and in L 2 (0, T ; X 1 0 (Ω)). Proposition 3.5. Assume that σ B 1 , σ B 2 , σ D , α, β, f 1 , f 2 , v 0 , s 0 are as in Subsection 2.3. Assume that I ion ≡ 0. Then, problem (2.26)-(2.34) admits a unique solution (V, U) ∈ L 2 (0, T ; H 1 null (Ω B )) × L 2 (0, T ; X 1 0 (Ω)), such that V ∈ C 0 ([0, T ]; L 2 (Ω B )) and [U] ∈ C 0 ([0, T ]; L 2 (Γ )).
As a consequence of the previous results, we finally get our main theorem.
Proof. The proof can be obtained following the same approach as in [13] (see, also, [14, §2.4.1]). Indeed, recalling that the function g appearing in the gating equation