Yves Bourgault

Yves Coudière
email: yves.coudiere@univ-nantes.fr

Charles Pierre
email: charles.pierre@univ-pau.fr

Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology ⋆

Keywords: reaction-diffusion equation, bidomain model, cardiac electrophysiology 1991 MSC: 35A05, 35K57, 35Q80

We study the well-posedness of the bidomain model, that is commonly used to simulate electrophysiological wave propagation in the heart. We base our analysis on a formulation of the bidomain model as a system of coupled parabolic and elliptic PDEs for two potentials and ODEs representing the ionic activity. We first reformulate the parabolic and elliptic PDEs into a single parabolic PDE by the introduction of a bidomain operator. We properly define and analyze this operator, basically a non differential and non local operator. We then present a proof of existence, uniqueness and regularity of a local solution in time through a semigroup approach, but that applies to fairly general ionic models. The bidomain model is next reformulated as a parabolic variational problem, through the introduction of a bidomain bilinear form. A proof of existence and uniqueness of a global solution in time is obtained using a compactness argument, this time for an ionic model reading as a single ODE but including polynomial nonlinearities. Finally, the hypothesis behind the existence of that global solution are verified for three commonly used ionic models, namely the FitzHugh-Nagumo, Aliev-Panfilov and MacCulloch models.

Introduction

The goal of the present paper is to investigate existence and uniqueness of solutions of the bidomain equations, commonly used for modeling the propagation of electrophysiological waves in the myocardium [START_REF] Colli Franzone | Spread of excitation in 3-D models of the anisotropic cardiac tissue. II: Effects of fiber architecture and ventricular geometry[END_REF][START_REF] Colli-Franzone | Modeling ventricular excitation: Axial and orthotropic anisotropy effects on wavefronts and potentials[END_REF][START_REF] Colli Franzone | Mathematical modeling of the excitation process in myocardial tissue: Influence of fiber rotation on wavefront propagation and potential field[END_REF][START_REF] Geselowitz | A bidomain model for anisotropic cardiac muscle[END_REF][START_REF] Gulrajani | A membrane-based computer heart model employing parallel processing[END_REF][START_REF] Henriquez | Simulating the electrical behavior of cardiac tissue using the bidomain model[END_REF][START_REF] Hooks | Cardiac microstructure: Implications for electrical propagation and defibrillation in the heart[END_REF][START_REF] Keener | A numerical method for the solution of the bidomain equations in cardiac tissue[END_REF][START_REF] Keener | Mathematical Physiology[END_REF][START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF][START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF][START_REF]Computational Biology of the Heart[END_REF][START_REF] Roth | Approximate analytical solutions to the bidomain equations with unequal anisotropy ratio[END_REF][START_REF] Roth | Mendering of spiral waves in anisotropic cardiac tissue[END_REF][START_REF] Vigmond | Computational tools for modeling electrical activity in cardiac tissue[END_REF][START_REF] Bourgault | Simulation of electrophysiological waves with an unstructured finite element method[END_REF]. The bidomain model was proposed thirty years ago [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF][START_REF] Schmidt | Biological information processing using the concept of interpenetrating domains[END_REF][START_REF] Tung | A bidomain model for describing ischemic myocardial D-C potentials[END_REF] but formal derivations of the model were obtained later [START_REF] Ambrosio | On the asymptotic behaviour of anisotropic energies arising in the cardiac bidomain model[END_REF][START_REF] Neu | Homogenization of syncitial tissues[END_REF][START_REF] Pennacchio | Multiscale modelling for the bioelectric activity of the heart[END_REF]. Consider a bounded subset Ω of R d (d = 2, 3) representing the myocadium. This model can be written as two degenerate parabolic PDEs with boundary conditions, coupled to a set of m ODEs, and some initial data:

∂u ∂t + f (u, w) -∇ • (σ i ∇u i ) = s i , in (0, +∞) × Ω, (1)
∂u ∂t + f (u, w) + ∇ • (σ e ∇u e ) = -s e , in (0, +∞) × Ω, (2)
∂w ∂t + g(u, w) = 0, in (0, +∞) × Ω, (3)
u = u i -u e , in (0, +∞) × Ω, (4) σ i ∇u i • n = 0, σ e ∇u e • n = 0, in (0, +∞) × ∂Ω,

(5) u(0) = u 0 , w(0) = w 0 , in Ω. [START_REF] Henriquez | Simulating the electrical behavior of cardiac tissue using the bidomain model[END_REF] Here, the unknowns are the functions u i (t, x) ∈ R, u e (t, x) ∈ R and w(t, x) ∈ R m (m ≥ 1), which are respectively the intra-and extra-cellular potentials and some ionic variables (currents, gating variables, concentrations, etc). The variable u defined in (4) denotes the transmembrane potential. Naturally, n denotes the unit normal to ∂Ω outward of Ω.

The other data, σ i,e (x) are conductivity matrices; f : R × R m → R and g : R×R m → R m are functions representing the ionic activity in the myocardium; and s i,e : (0, +∞) × Ω → R are external applied current sources. We point out that the conductivity matrices are quite specific: at each x ∈ Ω, they have the same eigenbasis, and for x ∈ ∂Ω, the normal n(x) to ∂Ω is an eigenvector of both σ i (x) and σ e (x). In an isolated heart, no current flows out of the heart, as expressed by the boundary conditions [START_REF] Gulrajani | A membrane-based computer heart model employing parallel processing[END_REF]. The initial data (at t = 0) is set only on u and w.

To our knowledge, only one proof of the well-posedness of the bidomain model is available in the literature [START_REF] Colli Franzone | Degenerate evolution systems modeling the cardiac electric field at micro and macroscopic level[END_REF]. This proof is based on a reformulation of ( 1)-( 6) as a Cauchy problem for an evolution variational inequality in a properly chosen Sobolev space. This approach apply only to the case f (u, w) = u(ua)(u -1) + w and g(u, w) = -ǫ(ku -w) known as the FitzHugh-Nagumo equations [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF]. Unfortunately, it is not well suited to modeling the action potential in myocardial excitable cells [9][12, chap. 1].

In this paper, we propose a new approach to the bidomain equations in order to address ionic models (i.e. functions f , g) adapted to the myocardial cells.

Our idea is to reformulate (1)-( 2) as a parabolic PDE coupled to an elliptic one, by replacing u i = u + u e in [START_REF] Colli Franzone | Spread of excitation in 3-D models of the anisotropic cardiac tissue. II: Effects of fiber architecture and ventricular geometry[END_REF] and substracting [START_REF] Colli-Franzone | Modeling ventricular excitation: Axial and orthotropic anisotropy effects on wavefronts and potentials[END_REF] to [START_REF] Colli Franzone | Spread of excitation in 3-D models of the anisotropic cardiac tissue. II: Effects of fiber architecture and ventricular geometry[END_REF]. The boundary condition ( 5) is also reformulated in terms of u and u e . The complete form of the new system is:

∂u ∂t + f (u, w) -∇ • (σ i ∇u) -∇ • (σ i ∇u e ) = s i , in (0, +∞) × Ω, (7)
∇ • (σ i ∇u + (σ i + σ e )∇u e ) = -(s i + s e ), in (0, +∞) × Ω, (8) ∂w ∂t + g(u, w) = 0, in (0, +∞) × Ω, ( 9)

σ i ∇u • n + σ i ∇u e • n = 0, in (0, +∞) × ∂Ω, (10)
σ i ∇u • n + (σ i + σ e )∇u e • n = 0, in (0, +∞) × ∂Ω, (11) u
(0) = u 0 , w(0) = w 0 , in Ω. ( 12
)
Hence we consider the problem of finding unknown functions u, u e and w verifying ( 7)- [START_REF]Computational Biology of the Heart[END_REF]. Naturally, the regular solutions of (1)-( 6) and ( 7)- [START_REF]Computational Biology of the Heart[END_REF] are exactly the same.

While the main difficulty in (1)-( 2) is the degeneracy in the temporal derivative, in our approach, we are able to replace u e in [START_REF] Hooks | Cardiac microstructure: Implications for electrical propagation and defibrillation in the heart[END_REF] by the solution of [START_REF] Keener | A numerical method for the solution of the bidomain equations in cardiac tissue[END_REF] with the boundary condition [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF], reducing ( 7)- [START_REF]Computational Biology of the Heart[END_REF] to a Cauchy problem for a single abstract parabolic equation with unknowns u and w, which reads

∂u ∂t + f (u, w) + Au = s, (13)
∂w ∂t + g(u, w) = 0, ( 14)

u(0) = u 0 , w(0) = w 0 , ( 15
)
where A is an integro-differential 2nd order elliptic operator accounting for the boundary conditions [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF] and [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF] and s is a modified source term, both given formally by

Au = -∇ • (σ i ∇u) + ∇ • σ i ∇ {∇ • (σ i + σ e )∇} -1 (∇ • σ i ∇u) , s = s i -∇ • σ i ∇ {∇ • (σ i + σ e )∇} -1 (s i + s e ) .
Afterwards, we are able to recover u e and u i with

u i = u + u e , u e = {∇ • (σ i + σ e )∇} -1 (-(s i + s e ) -∇ • σ i ∇u) . (16)
The unknowns u i and u e are defined up to an additional constant.

Two specific definitions of the so-called bidomain operator A and source term s will be given in §3, definitions 5 and 12, in order to express ( 13)- [START_REF] Vigmond | Computational tools for modeling electrical activity in cardiac tissue[END_REF] in strong and weak variational senses.

We choose to study both strong (see §4, definition 12) and weak (see §5, definition 5 and lemma 9) solutions of ( 7)-( 12), using the operator A and the source term s. The strong solution theory applies to fairly general ionic models [START_REF] Beeler | Reconstruction of the Action Potential of Ventricular Myocardial Fibres[END_REF][START_REF] Di-Francesco | Simulations of Ionic Currents and Concentration Changes[END_REF][START_REF] Luo | A model of the Ventricular Cardiac Action Potential[END_REF][START_REF] Noble | A modification of the Hodgkin-Huxley equation applicable to Purkinje fibre action and pacemaker potentials[END_REF], but only provides solutions that are local in time. On the other hand, weak solutions are obtained for simpler ionic models reading as a single ODE with polynomial nonlinearities, but are global in time. These ionic models include the FitzHugh-Nagumo model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF] and simple models more adapted to myocardial cells, such as the Aliev-Panfilov [START_REF] Panfilov | A simple two-variable model of cardiac excitation[END_REF] and MacCulloch [START_REF] Roger | A collocation-Galerkin finite element model of cardiac action potential propagation[END_REF] models.

Strong solutions on an interval [0, τ ) are functions u(t), u e (t), w(t) with value respectively in H 2 (Ω), H 2 (Ω)/R, L ∞ (Ω) such that ( 7)-( 9) holds for all t ∈ (0, τ ) and a.e. x ∈ Ω and ( 10)-( 11) holds for all t ∈ (0, τ ) and a.e. x ∈ ∂Ω (definition 18). There is also a notion of mild solutions, but these coincide with strong solutions because we are interested in continuous solutions. See [START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF] for details.

As correctly stated in definition 26, weak solutions on an interval [0, τ ) are functions u(t), u e (t), w(t) with value respectively in

H 1 (Ω), H 1 (Ω)/R, L 2 (Ω) such that d dt (u(t), v) + Ω σ i ∇(u(t) + u e (t)) • ∇v + Ω f (u(t), w(t))v = Ω s i (t)v d dt (w(t), v) + Ω g(u(t), w(t))v = 0
respectively for all v ∈ H 1 (Ω) and v ∈ L 2 (Ω), where for a.e. t ∈ (0, τ ),

Ω σ i ∇u(t) • ∇v e + Ω (σ i + σ e )∇u e (t) • ∇v e = Ω (s i (t) + s e (t))v e
for all v e ∈ H 1 (Ω)/R. They are equivalently defined as weak solutions to ( 13)-( 14) with u e given by ( 16) in a weak sense (lemma 28).

Our main results are :

• Existence, uniqueness and regularity for strong solutions under weak assumptions on f and g, but only with τ > 0 small enough, and for regular data (∂Ω, σ i,e , s i,e and u 0 , w 0 );

• Existence for weak solutions under more restrictive assumptions on f and g, but for τ = +∞ and with minimal regularity on the data.

The first important difficulty concerns the definition of the mappings (t, x) → f (u(t, x), w(t, x)) and (t, x) → g(u(t, x), w(t, x)):

• Either f : R × R m → R, g : R × R m → R m are locally Lipschitz functions and u, w are regular enough and we can find τ > 0 small enough for strong solutions to be defined on [0, τ ); • Or f, g are only continuous but with polynomial growth at infinity with suitable power (as given by assumption (H3), §5) which allow for weak solutions to exist for any t > 0.

The second difficulty concerns the operator A. Clearly, it is uniformly elliptic. In the semigroup approach, the solutions are obtained through a fixed point technique, but only with τ > 0 small enough. Since the functions f, g of interest usually define vector fields with invariant regions, as defined in [START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF] (think of the FitzHugh-Nagumo equations), global solutions (τ = +∞) are expected from a maximum principle on A. It is interesting to note that global existence results for system of reaction-diffusion equations, similar to our bidomain model, are only available when the elliptic operators in the system is a Laplacian or more generally a second order elliptic operator reading the form ∇ • (σ∇•) [START_REF] Amour | The Cauchy problem for a coupled semilinear parabolic system[END_REF][START_REF] Dickstein | A maximum principle for semilinear parabolic systems and applications[END_REF][START_REF] Guarguaglini | Global existence and uniqueness of solutions for multidimensional weakly parabolic systems arising in chemistry and biology[END_REF][START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF][START_REF] Petersson | On global existence for semilinear parabolic systems[END_REF][START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF][START_REF] Snoussi | Global existence, asymptotic behavior and self-similar solutions for a class of semilinear parabolic systems[END_REF][START_REF] Temam | Infinite-Dimensional Dynamical Systems in Mechanics and Physics[END_REF]. In these cases, essential properties of elliptic second order differential operators are required to derive comparison theorems or maximum principles. Unfortunately, we were not able to prove any comparison theorem, because of the integro-differential form of our operator A. And our strong solutions are proved to exist only for t < τ , with τ small enough.

There are numerical experiments and intuitions that the solution are bounded functions, hence solutions should exist for all times. Here, the variational approach is helpful to obtain solutions for any t > 0 (τ = +∞). Of course, it requires an energy estimates, that depends on an additional assumption on the function f, g (besides their polynomial growth) stated as (H4). The variational approach is quite technical, for sake of simplicity it is developped only for m = 1.

The paper is sketched as follows: some notations and explanations on the boundary conditions are introduced in §2; the bidomain operator is defined in §3; results on local in time strong solutions are given in §4; the variational approach is covered in §5, and the three examples are considered in §6.

In §4, we shall briefly but completely define and prove existence, uniqueness and regularity results for strong solutions on a local time interval. The analysis mainly explain how the problem enters the framework of solutions as defined in [START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF]. The results are obtained for smooth data.

In §5, our main result is stated: existence of a weak solutions. It is based on a classical Galerkin technique, a priori estimates and compactness results [START_REF] Lions | Résolution de quelques problèmes aux limites non-linéaires[END_REF].

The results are obtained under minimal regularity assumptions on the data. The problem of regularity of the solutions is not addressed here.

Notations and Boundary Conditions

In order to formulate the weak form of the equations, we only need Ω to be a bounded open subset of R d with Lipschitz boundary ∂Ω, and the conductivity matrices σ i,e to be functions of the space variable x ∈ Ω with coefficients in L ∞ (Ω) and uniformly elliptic. Namely, we assume that there exists constants 0 < m < M such that

m |ξ| 2 ≤ ξ t σ i,e (x) ξ ≤ M |ξ| 2 , ∀ξ ∈ R d , (17)
for a.e. x ∈ Ω. The matrices σ i,e are symmetric. Now, to define strong solutions and for the boundary conditions to make sense, we need additional regularity on these data: ∂Ω is supposed to be C 2+ν and σ i,e to have their coefficients in C 1+ν ( Ω) for some ν ≥ 0. Of course condition ( 17) is supposed to be valid for any x ∈ Ω.

The source terms s i (t) and s e (t) are related to the applied stimulating currents. Their regularity will be specified in §4 and §5. We point out that s i + s e must have a zero mean value. The physical reason for this is that there is no current flowing out of the heart through its boundary as stated by boundary conditions [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF] or [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF], and that the intra-and extra-cellular medias are electrically communicating through the cells membrane. Therefore the current conservation applied to the whole heart, in other words Eq. ( 8) integrated over all Ω ensures that s i + s e has a zero mean value. That is:

Ω (s i (x) + s e (x)) dx = 0. ( 18
)
We now emphasize on a specificity of the bidomain equations:

(1) the symmetric matrices σ i,e (x) have the same basis of eigenvectors Q(x) = (q 1 (x), . . . , q d (x)) in R d , which reflect the organization of the muscle in fibers (direction q 1 (x)) [START_REF] Colli Franzone | Mathematical modeling of the excitation process in myocardial tissue: Influence of fiber rotation on wavefront propagation and potential field[END_REF]; and then σ i,e (x) = Q(x)Λ i,e (x)Q(x) T where Λ i,e (x) = diag(λ 1 i,e (x), . . . , λ d i,e (x)). (2) the muscle fibers are tangent to the boundary ∂Ω [START_REF] Colli Franzone | Mathematical modeling of the excitation process in myocardial tissue: Influence of fiber rotation on wavefront propagation and potential field[END_REF], so that σ i,e (x)n(x) = λ d i,e n(x), a.e. x ∈ ∂Ω [START_REF] Tung | A bidomain model for describing ischemic myocardial D-C potentials[END_REF] with λ d i,e (x) ≥ m > 0.

As a consequence, we have the Lemma 1 Suppose that Ω has a C 1 boundary ∂Ω, Q(x) and Λ i,e (x) are C 0 ( Ω).

For functions u, u e ∈ H 2 (Ω), the conditions (5), and the conditions (10)- [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF], and the homogeneous Neumann conditions

∇u i • n = 0, ∇u e • n = 0, in ∂Ω are equivalent.
Proof. Conditions ( 5) and ( 10)-( 11) are linear combination one of the other, equivalent to the conditions above because of identity [START_REF] Tung | A bidomain model for describing ischemic myocardial D-C potentials[END_REF].

The bidomain operator

Now we need to study the mapping

(u, s i + s e ) → solve (8)
u e → replace in ( 7)

s i + ∇ • (σ i ∇u) + ∇ • (σ i ∇u e ) := Au + s.
Hence, we study the system

-∇ • (σ i ∇u) -∇ • (σ i ∇u e ) = s i , in Ω, ( 20
) ∇ • (σ i ∇u + (σ i + σ e )∇u e ) = -(s i + s e ), in Ω, ( 21
)
with the boundary conditions [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF] and [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF], given conductivity matrices σ i (x), σ e (x) and data s i (x), s e (x) such that Ω (s i + s e ) = 0 (or more generally s i + s e , 1 = 0).

Variational formulation

Equations ( 20) and ( 21) with the Homogeneous Neumann boundary conditions [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF] and [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF] have solutions (u, u e ) defined only up to an additive constant.

In practice, the nonlinear term determines u in H 1 (Ω) but u e remains defined up to an additive constant. Weak solutions will be found in H 1 (Ω)/R, using a bilinear form in

H 1 (Ω)/R × H 1 (Ω)/R, that is extended to H 1 (Ω) × H 1 (Ω)
in order to address the bidomain equations.

Given a Banach space X of functions integrable on Ω, we define its subspace X/R = {u ∈ X, Ω u = 0} ⊂ X, that is a Banach space with the norm u X/R = u X . For any u ∈ X, we also note

[u] = u -Ω u ∈ X/R.
In practice, we are working with V = H 1 (Ω), H = L 2 (Ω) endowed with their usual norms. We shall note U = V /R, and we have

U ⊂ H/R ≡ (H/R) ′ ⊂ U ′ , D(Ω) ⊂ V ⊂ H ≡ H ′ ⊂ V ′ ⊂ D ′ (Ω)
where the injections are continuous and the injection V → H and U → H/R are compact. Note that |u| U = ( Ω |∇u| 2 )

1/2 defines a norm on U equivalent to the norm on V , due to the Poincaré-Wirtinger inequality

∃C > 0, ∀u ∈ U, Ω |u| 2 ≤ C|u| 2 U := C Ω |∇u| 2 . ( 22
)
Hence we shall use

| • | U instead of the norm induced by • V in U . An element s in H/R or H is identified to the element s ∈ (H/R) ′ by s, v = Ω sv for all v ∈ H/R. An element s ∈ V ′ is
identified to an element of U ′ by just restricting the duality product s, v := s, v V ′ ×V to functions v in the subspace U of V . Conversely, an element s ∈ U ′ can only be extended to the whole space V independently of the value v -[v] using a special condition like s, 1 = 0.

Given a regular solution u, u e ∈ H 2 (Ω)/R, we obtain a variational formulation by multiplying Eq. ( 20) by a test function v ∈ U , multiplying Eq. ( 21) by a test function v e ∈ U , integrating by parts the second order terms and adding the two resulting variational equations. The boundary terms vanish due to the boundary conditions ( 10) and ( 11) on u and u e . The resulting variational problem reads: Find (u, u e ) ∈ U × U such that

a i (u, v) + a i (u e , v) + a i (u, v e ) + (a i + a e )(u e , v e ) = s i , v + s i + s e , v e , (23)
for all (v, v e ) ∈ U × U . The bilinear forms a i,e (•, •) on U × U are

a i (u, v) = Ω σ i ∇u • ∇v dx, a e (u, v) = Ω σ e ∇u • ∇v dx, ∀(u, v) ∈ U × U,
and s i , (s i + s e ) ∈ V ′ are given source terms.

Under the hypothesis (17), the bilinear forms a i,e (•, •) are symmetric continuous and uniformly elliptic on U . We define the bilinear form

b((u, u e ), (v, v e )) = a i (u, v) + a i (u e , v) + a i (u, v e ) + (a i + a e )(u e , v e )
on U × U , for the sake of simplicity. We have

Lemma 2 The bilinear form b(•, •) is symmetric, continuous and uniformly elliptic on (U × U ) × (U × U ) for the norm |(v, v e )| U ×U = max(|v| U |, |v e | U ).
Proof. Clearly, the bilinear form b(•, •) can be rewritten as

b((u, u e ), (v, v e )) = a i (u + u e , v + v e ) + a e (u e , v e ).
It is bilinear and symmetric. Additionally, a simple computation shows that

1 3 |(u, u e )| 2 U ×U ≤ |u + u e | 2 U + |u e | 2 U ≤ 6|(u, u e )| 2
U ×U , and using the inequalities (17) we have

|b((u, u e ), (v, v e ))| ≤ 6M |(u, u e )| U ×U |(v, v e )| U ×U , b((u, u e ), (u, u e )) ≥ m 3 |(u, u e )| 2 U ×U ,
which proves that b is continuous and coercive.

Theorem 3 Let s i , s e ∈ V ′ and u ∈ U be given. The variational equations

(a i + a e )(ũ e , v e ) = -a i (u, v e ), ∀v e ∈ U, (24)
(a i + a e )(ū e , v e ) = s i + s e , v e , ∀v e ∈ U, (25)
have unique solutions ũe , ūe ∈ U . For any u, v ∈ U , we can define the mappings

ā(u, v) = b((u, ũe ), (v, 0)), s, v = s i , v -a i (ū e , v). ( 26
)
The mapping ā is bilinear, symmetric, continuous and uniformly elliptic on U × U , and the mapping s is linear and continuous on U .

Eq. ( 23) has a unique solution (u, u e ) where u is also the unique solution of

ā(u, v) = s, v , ∀v ∈ U, (27)
and u e = ũe + ūe , where ũe , ūe are the solutions of ( 24) and (25).

Proof. The idea is that, for any (u,

u e ) ∈ U × U , b((u, u e ), (v, v e )) = s i , v + s i + s e , v e , ∀(v, v e ) ∈ U × U, (28)
⇔ b((u, u e ), (v, 0)) = s i , v , ∀v ∈ U, (29) b
((u, u e ), (0, v e )) = s i + s e , v e , ∀v e ∈ U, (30)
⇔ b((u, ũe ), (v, 0)) = s i , v -b((0, ūe ), (v, 0)), ∀v ∈ U, (31) b
((u, ũe ), (0, v e )) = 0, ∀v e ∈ U, (32) b
((0, ūe ), (0, v e )) = s i + s e , v e , ∀v e ∈ U, (33)
u e = ũe + ūe . (34)
Eqs. ( 32) and ( 33) are exactly Eqs. ( 24) and ( 25), respectively; and Eq. ( 31) is exactly [START_REF] Luo | A model of the Ventricular Cardiac Action Potential[END_REF] with the notations [START_REF] Di-Francesco | Simulations of Ionic Currents and Concentration Changes[END_REF].

Note that the mappings v e ∈ U → a i (u, v e ) ∈ R and v e ∈ U → s i +s e , v e ∈ R are linear and continuous. Equations ( 24) and ( 25) have unique solutions ũe and ūe by the theorem of Lax-Milgram. We have

m 2M |u| U ≤ |ũ e | U ≤ M 2m |u| U , |ū e | U ≤ C 2m s i + s e U ′ , (35)
where C > 0 is the constant of the Poincaré-Wirtinger inequality [START_REF] Pennacchio | Multiscale modelling for the bioelectric activity of the heart[END_REF]. The mapping s defined in ( 26) is obviously linear, and continuous because:

∀v ∈ U, | s, v | ≤ C s i U ′ + CM 2m s i + s e U ′ |v| U .
It remains to prove that the mapping ā defined in ( 26) is bilinear on U × U , continuous, uniformly elliptic and symmetric. Consider the function ṽe ∈ U constructed like ũe (the solution of ( 32)) as the solution of b((v, ṽe ), (0, v e )) = 0 for all v e ∈ U . The result is deduced from the definition of b(•, •), lemma 2 and (35) because

ā(u, v) = b((u, ũe ), (v, 0)) = b((u, ũe ), (v, ṽe )) = b((v, ṽe ), (u, ũe )) = b((v, ṽe ), (u, 0)) = ā(v, u).
Specifically, we have

|ā(u, v)| ≤ M 1 + M 2m |u| U |v| U , ā(u, u) ≥ m 3 1 + m 2M |u| 2 U .
Remark 4 Conversely, let u be the solution to [START_REF] Luo | A model of the Ventricular Cardiac Action Potential[END_REF] and u e = ũe + ūe be given by ( 24) and (25). The space U does not contain the space D(Ω), but for any v ∈ D(Ω), we have ∇v = ∇[v] and [v] ∈ U . We must impose the extra condition s i + s e , 1 = 0 to get

s i + s e , v = s i + s e , [v] .
In that case u e = ũe + ūe can easily be proved to verify

∇ • (σ i ∇u + (σ i + σ e )∇u e ) = -(s i + s e ), in D ′ (Ω).
Additionally, if u, u e ∈ H 2 (Ω), then the equation is verified a.e. in Ω and the boundary condition (11) holds a.e. in ∂Ω At last, the operator ā is extended to V ×V , in order to state the full bidomain problem.

Definition 5

The bidomain bilinear form a is defined on

V × V by a(u, v) = ā([u], [v]), ∀(u, v) ∈ V × V.
Given s i , s e ∈ V ′ such that s i + s e , 1 = 0, the source term s is extended to a linear form on V , still denoted by s, by

s, v = s i , v -a i (ū e , [v]), ∀v ∈ V,
where ūe is given by [START_REF] Beeler | Reconstruction of the Action Potential of Ventricular Myocardial Fibres[END_REF].

Theorem 6

The bilinear form a(•, •) is symmetric, continuous and coercive on V ,

∀u ∈ V, α u 2 V ≤ a(u, u) + α u 2 H , (36) ∀u, v ∈ V, |a(u, v)| ≤ M u V v V , (37)
with α = m 3 1 + m 2M and M = 3M 1 + M 2m .
There exists an increasing sequence 0 = λ 0 < . . . ≤ λ i ≤ . . . in R and an orthonormal Hilbert basis of H of eigenvectors (ψ i ) i∈N such that for all i ∈ N,

ψ i ∈ V and ∀v ∈ V, a(ψ i , v) = λ i (ψ i , v). ( 38
)
Given s i , s e ∈ V ′ such that s i , 1 = s e , 1 = 0, if u ∈ V is a solution to a(u, v) = s, v , ∀v ∈ V, (39)
and u e = ũe + ūe ∈ U is given by ( 24) and (25), then (u, u e ) is a weak solution to (20), ( 21) with the boundary conditions (10), [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF]. If additionally u, u e ∈ H 2 (Ω), then they verify (20), ( 21) a.e. in Ω and (10), ( 11) a.e. in ∂Ω.

Proof. The bilinear form a(•, •) is well-defined and symmetric. It verifies [START_REF] Petersson | On global existence for semilinear parabolic systems[END_REF] and ( 37) because of theorem 3, and because

u V = ( u 2 H + |[u]| 2 U ) 1/2 .
The spectral problem is very classical, and can be found for instance in [40, th. 6.2-1 and rem. 6.2-2, p.137-138]. In this case, we have λ 0 = 0 because the bilinear form a obviously vanishes only for constant functions.

The equivalence with the strong formulation and the boundary conditions is also a classical result, stated partly in remark 4, and easily deduced from the definition 5 of a and s:

a(u, v) = s, v ⇔    b((u, u e ), ([v], 0)) = s i , v ∀v ∈ V, b((u, u e ), (0, v e )) = s i + s e , v e ∀v e ∈ U.
The second equation is equivalent to [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF] and [START_REF] Beeler | Reconstruction of the Action Potential of Ventricular Myocardial Fibres[END_REF] to define u e from [u] ∈ U and the first equation is obviously the weak form of [START_REF] Ambrosio | On the asymptotic behaviour of anisotropic energies arising in the cardiac bidomain model[END_REF].

Remark 7 The two conditions s i , 1 = s e , 1 = 0 are needed for the solution u ∈ U = V /R of (39) to be interpreted as a weak solution in V of the partial differential equations (20), [START_REF] Neu | Homogenization of syncitial tissues[END_REF] with the Homogeneous Neuman boundary conditions [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF] and [START_REF] Lines | Modeling the electrical activity of the heart: A bidomain model of the ventricles embedded in a torso[END_REF].

For the full nonlinear bidomain problem, only s i + s e , 1 = 0 will be required with no zero-average condition on s i and s e taken individually. The definition 5 requires only that s i + s e , 1 = 0.

Weak operator

This variational process can also be handled through operators. These are weak operators, defined from U onto U ′ or V onto V ′ . Indeed we are able to define A i,e and A by duality by setting

A i,e u, v = a i,e (u, v), Au, v = ā(u, v), ∀(u, v) ∈ U × U.
They are all one-to-one continuous mappings from U onto U ′ , with continuous inverse (from Lax-Milgram theorem).

Lemma 8 Given s i , s e ∈ U ′ , the source term s ∈ U ′ defined in theorem 3 is such that

s = s i -A i (A i + A e ) -1 (s i + s e ) = -s e + A e (A i + A e ) -1 (s i + s e ),
and we have

A = A i (A i + A e ) -1 A e = (A -1 i + A -1 e ) -1 .
Proof. We can rewrite ( 24) and ( 25) as

(A i + A e )ū e = s i + s e , (A i + A e )ũ e = -A i u.
The result is found by substituting ũe , ūe in [START_REF] Di-Francesco | Simulations of Ionic Currents and Concentration Changes[END_REF]. A short computation shows that

A = A i -A i (A i + A e ) -1 A i = A i (Id -(A i + A e ) -1 A i ) = A i (A i + A e ) -1 (A i + A e -A i ).
The second equality defining A comes from a simple algebraic manipulation and the fact that the operator

A -1 i + A -1 e is invertible.
The second formulation of s is due to the identity

s -s = s i + s e -(A i + A e )(A i + A e ) -1 (s i + s e ) = 0. Lemma 9 Define J : u ∈ V → [u] ∈ U and its transpose J T : U ′ → V ′ .
For any s i , s e ∈ V ′ with s i + s e , 1 = 0, the source term s ∈ V ′ and the bilinear operator a given by definition 5 are such that

s = s i -J T A i (A i + A e ) -1 (s i + s e ) = -s e + J T A e (A i + A e ) -1 (s i + s e ), and
a(u, v) = Au, v for all (u, v) ∈ V × V, with A = J T AJ : V → V ′ . Proof. It is obvious since a(u, v) = ā([u], [v]) = AJu, Jv = J T AJu, v for all (u, v) ∈ V × V . Concerning s, we have for all v ∈ V , s, v = s i , v -a i (ū e , [v]) = s i , v -A i (A i + A e ) -1 (s i + s e ), Jv = s i -J T A i (A i + A e ) -1 (s i + s e ), v .
Remark 10 Although it is a positive operator, A is not in general a differential operator, being the harmonic average of two second order differential elliptic operators.

f ∈ H/R we have u = (A i,e ) -1 f ∈ H 2 (Ω)/R, or u = (A i + A e ) -1 f ∈ H 2 (Ω)/R.
As a consequence, we have

∇ • (σ i,e ∇u) = f a.e. in Ω, σ i,e ∇u • n = 0 a.e. in ∂Ω, ∇ • ((σ i + σ e )∇u) = f a.e. in Ω, (σ i + σ e )∇u • n = 0 a.e. in ∂Ω.
With lemma 1 in addition, we can define unbounded operators in H/R, still denoted by A i and A e and A i + A e , with domains [START_REF] Raviart | Introduction à l'analyse numériaue des équations aux dérivées partielles[END_REF] with

D(A i ) = D(A e ) = D(A i + A e ) = D(A)/R by A i u = ∇•(σ i ∇u), A e u = ∇•(σ e ∇u), (A i +A e )u = ∇•((σ i +σ e )∇u),
D(A) = u ∈ H 2 (Ω), ∇u • n = 0 a.e. in ∂Ω ⊂ H. (41)

Lemma 11

The operators A i , A e , A i + A e are maximal monotone in H/R and self-adjoint. They have compact inverses in H/R.

Proof. The operators

A i , A e , A i +A e verify, for all (u, v) ∈ D(A)/R×D(A)/R, A i,e u, v = a i,e (u, v), (A i + A e )u, v = (a i + a e )(u, v).
They are obviously maximal-monotone and self-adjoint. Their inverses

(A i,e ) -1 : H/R → H/R, (A i + A e ) -1 : H/R → H/R
are compact operators since their range is D(A)/R and the injection D(A)/R → H/R is compact.

Definition 12

Given s i,e ∈ H such that s i + s e ∈ H/R, we define the strong bidomain operator A : D(A) ⊂ H → H and the source term s ∈ H by:

Au = A i (A i + A e ) -1 A e [u], ∀u ∈ D(A), (42) s
= s i -A i (A i + A e ) -1 (s i + s e ) = -s e + A e (A i + A e ) -1 (s i + s e ). ( 43
)
The definitions make sense because [u] ∈ D(A)/R = D(A i ) and the range of

(A i + A e ) -1 is D(A)/R = D(A e )
, and of course H/R ⊂ H. We use the same notation for the weak and strong bidomain operators as it will be clear from context which operator we refer to.

Theorem 13 Consider s i , s e ∈ H such that s i + s e ∈ H/R. The strong bidomain operator A of definition 12 is self-adjoint and maximal-monotone. We have

∀(u, v) ∈ D(A) × D(A), (Au, v) = a(u, v),
and the source term s ∈ H of definition 12 can be identified to the source term s ∈ V ′ of definition 5 through the identity s, v = (s, v) for all v ∈ V ⊂ H.

The sequence (λ i ) i∈N and the orthonormal Hilbert basis (ψ i ) i∈N of eigenvectors defined in theorem 6 is such that ψ i ∈ D(A) for all i ∈ N and

D(A) = {u ∈ H, i≥0 λ 2 i (u, ψ i ) 2 < ∞}, Au = i≥0 λ i (u, ψ i )ψ i . (44)
For u ∈ H, we have Au = s and u e = (A i + A e ) -1 ((s i + s e ) -A i u) ∈ D(A)/R ⇔ (u, u e ) verify (20) and (21) a.e. in Ω and the boundary conditions [START_REF] Lines | Mathematical models and numerical methods for the forward problem in cardiac electrophysiology[END_REF] and

(11) a.e. in ∂Ω.

Proof. For u ∈ D(A), Au ∈ H/R ⊂ H is well defined. Consider ũe = -(A i + A e ) -1 A i [u] ∈ D(A)/R and ūe = (A i + A e ) -1 (s i + s e ) ∈ D(A)/R .
They are exactly the solutions to [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF] and [START_REF] Beeler | Reconstruction of the Action Potential of Ventricular Myocardial Fibres[END_REF]. We have

Au = A i [u]-A i (A i +A e ) -1 A i [u]
(simple computation) for all u ∈ D(A) and then for all v ∈ D(A),

(Au, v) = (A i [u], [v]) + (A i ũe , [v]) = a i ([u], [v]) + a i (ũ e , [v]) = b(([u], ũe ), (v, 0)) = a(u, v), (s, v) = (s i , v) -(A i ūe , [v]) = (s i , v) -a i (ū e , [v]) = s, v .
The remaining results are deduced from theorem 6. A is maximal-monotone (coercivity of a) and self-adjoint (symmetry of a) in H. The eigenvectors

ψ i ∈ V = H 1 (Ω) are such that Aψ i = λ i ψ i in V ′ and then ψ i ∈ D(A) (regularity result). As a consequence Aψ i = λ i ψ i in H and (44) is valid. The equivalence is true because if u ∈ D(A) ⊂ H 2 (Ω) then Au = s ⇔ a(u, v) = (s, v) ∀v ∈ V
and u e = (A i + A e ) -1 ((s i + s e ) -A i u) = ũe + ūe where ũe and ūe are the solutions to [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF] and [START_REF] Beeler | Reconstruction of the Action Potential of Ventricular Myocardial Fibres[END_REF].

Strong Solutions

The existence and uniqueness of strong solutions for ( 13)-( 15) is established in the framework of analytical semigroups, as presented in D. Henry's monograph [START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF].

Specific Assumptions and Notations

The result for the existence and uniqueness of strong solutions holds if we assume that Ω has a C2 boundary ∂Ω and that σ i,e have C 1 ( Ω) coefficients, in order to apply the definition and lemma from §3.3.

The integer m can be chosen as big as one wishes. As precisely stated below, the reaction terms f and g are only assumed to be locally Lipschitz functions (for existence and uniqueness of solutions results), so covering a wide range of realistic ionic models [START_REF] Beeler | Reconstruction of the Action Potential of Ventricular Myocardial Fibres[END_REF][START_REF] Di-Francesco | Simulations of Ionic Currents and Concentration Changes[END_REF][START_REF] Luo | A model of the Ventricular Cardiac Action Potential[END_REF][START_REF] Noble | A modification of the Hodgkin-Huxley equation applicable to Purkinje fibre action and pacemaker potentials[END_REF].

The second unknown w (a m-dimensional vector here) will be searched in a Banach space B m = B × B . . . × B where

• either B = L ∞ (Ω),
• or B = C ν (Ω), the space of all (globally) ν-Hölder continuous functions on Ω 2 . This last choice will be needed to establish the regularity of the solutions. In the sequel, ν will represent a real number 0 < ν < 1.

The function f : R × R m → R and g : R × R m → R m are:

• locally Lipschitz continuous functions on R × R m when assuming that B = L ∞ (Ω), • C 2 (R × R m ) regular functions when assuming that B = C ν (Ω).
At last, we assume the functions s i , s e to be locally ν-Hölder continuous in time, s i,e ∈ C ν loc ([0, +∞), H) for some ν > 0:

∀[t 1 , t 2 ] ⊂ [0, +∞), ∃C > 0, δ 1 , δ 2 ∈ [t 1 , t 2 ] ⇒ s i,e (δ 1 ) -s i,e (δ 2 ) H ≤ C|δ 1 -δ 2 | ν . ( 45
)
Consider Z = H × B m , with the norm (u, w) Z = max( u H , w B m ). It is a Banach space. We introduce the unbounded operator A in Z defined by

A : D(A) ⊂ Z → Z, Az = (Au, 0) ∈ Z, for z = (u, w) ∈ Z, with D(A) = D(A)
s(δ 1 ) -s(δ 2 ) = s i (δ 1 ) -s i (δ 2 ) -A i (A i + A e ) -1 (s i (δ 1 ) -s i (δ 2 ) + s e (δ 1 ) -s e (δ 2 )).
The result is obvious since

A i (A i + A e ) -1 : D(A)/R → D(A)/R is bounded.
Our next problem is to define the mapping

F : (u, w) ∈ Z → (f (u, w), g(u, w)) ∈ Z.
To get rid of that difficulty, we introduce the fractional powers A α and the interpolation spaces Z α = D(A α ). For α ≥ 0 the unbounded operator A α : D(A α ) ⊂ Z → Z is defined by:

Z α = {u ∈ H, i≥0 λ 2α i (u, ψ i ) 2 < ∞} × B m , A α (u, w) =   i≥0 λ α i (u, ψ i )ψ i , 0   .
The spaces Z α , equipped with the norm u α = u + A α u Z , are Banach spaces. Moreover (see [START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF]), for any 0 ≤ α ≤ β, we have the continuous and dense embedding Z β ⊂ Z α . These spaces form a sequence of decreasing functional spaces composed of functions whose regularity increases from Z

(α = 0) to D(A) ⊂ H 2 (Ω) × B m (α = 1)
. With the regularity we assumed for ∂Ω, we have the following embedding lemma:

Lemma 16 (Case B = L ∞ (Ω)) For B = L ∞ (Ω), f, g locally Lipschitz con- (4'
) for all t ∈ (0, τ ), u(t) ∈ D(A), u e (t) ∈ D(A)/R, and z verify for all t ∈ (0, τ ),

dz dt (t) + Az(t) + F (z(t)) = 0 in Z, (46)
using the previous definitions of A and F , while u e is given by

u e (t) = (A i + A e ) -1 ((s i (t) + s e (t)) -A i [u(t)]) ∈ D(A)/R. ( 47
)
Theorem 20 (Local existence and uniqueness) Consider 0 < α < 1 defined by lemma 16 (case B = L ∞ ) or lemma 17 (case B = C ν ) such that F : Z α → Z is well-defined and locally Lipschitzian. Then for any (u 0 , w 0 ) ∈ Z α , there exists T > 0 such that the problem (7)-( 12) has a unique strong solution on [0, T ) in the sense of definition 18.

We point out that choosing α such that F : Z α → Z is locally Lipschitzian imposes a strong constraint on the initial data z 0 = (u 0 , w 0 ) ∈ D(A α ). In dimension d = 3 for instance, one must have α > 3/4.

Proof. That theorem is a direct application of the local existence and uniqueness theorem in [31, p.54] since:

• there always exists 0 ≤ α < 1 such that F extend to a function F : Z α → Z locally Lipschitzian, for d = 1, 2, 3. • A is sectorial (lemma 14), • t → S(t) is locally ν-Hölder continuous for some ν > 0.

Regularity of the solutions

Given a real number 0 < ν < 1, we will assume throughout that subsection that B = C ν (Ω) and that the reaction terms f and g have C 2 regularity on R × R m . We will moreover assume that the boundary ∂Ω of the domain has C 2+ν regularity, and that σ i,e have their coefficients in C 1+ν ( Ω).

The operator A has a smoothing effect on the solutions of (46): for an initial data u 0 ∈ D(A α ), the solution satisfy u(t) ∈ D(A) for t > 0. This is due to the following elliptic regularity result (see [42, p.128] or [41, p.182]):

Lemma 21 Let σ be a uniformly elliptic tensor on Ω whose components belong to C 1+ν ( Ω) for some ν > 0. We also assume the boundary ∂Ω to have

C 2+ν regularity. If u ∈ D(A) satisfies ∇ • (σ∇u) ∈ C ν (Ω), then u ∈ C 2+ν (Ω).
Beyond this space smoothing effect on the first unknown u, some regularity in time also takes place [31, p.71

]:

Lemma 22 Let z : t ∈ (0, T ) → z(t) ∈ D(A) = Z 1 be the solution of the Cauchy problem (46) given by theorem 20. We have Z 1 ⊂ Z ν for any ν ≤ 1, and the solution moreover satisfies: t ∈ (0, T ) → z(t) ∈ Z ν is continuously (Fréchet) differentiable for any ν < 1.

Together with the previous elliptic regularity argument, these two results imply that the solutions for (46) actually are classical solutions provided that the initial data w 0 for the second variable w is smooth enough.

Theorem 23 (Regularity of the strong solution) Consider d/4 < α < 1 and 0 < ν < 2α -d/2, and assume that s i,e : [0, +∞) → H are locally ν-Hölder continuous and such that s i,e (t) ∈ C ν (Ω) for all t ≥ 0. For z 0 = (u 0 , w 0 ) ∈ Z α the unique solution z of (46) defined on [0, T ) for some T > 0 satisfies furthermore:

• Given any x ∈ Ω, the function t ∈ (0, T ) → z(t, x) is continuously differen- tiable in t, • Given any t ∈ (0, T ), the function x ∈ Ω → u(x, t) is twice continuously differentiable in x, i.e. u(•, t) ∈ C 2 (Ω).
Proof. Using the embedding from lemma 17 ensures that the solution t ∈ (0,

T ) → z(t) ∈ C ν (Ω) × (C ν (Ω)) m is continuously (Fréchet) differentiable.
This actually implies that (t, x) ∈ (0, T ) × Ω → z(x, t) = (u(x, t), w(x, t)) is continuously differentiable in t.

Let us now prove that Au(t) ∈ C ν (Ω) for t ∈ (0, T ). One has Au(t) = -du/dt(t) -f (u(t), w(t)) + s(t). Easily, we have f (u(t), w(t)) ∈ C ν (Ω) and also du/dt(t) ∈ C ν (Ω) thanks to lemma 22. Now s(t) = -s e (t) + A e (A i + A e ) -1 (s i (t) + s e (t)) and (s i + s e )(t) ∈ C ν (Ω) by assumption. By lemma 21, the function (A i + A e ) -1 (s i (t + s e (t)) belongs to C 2+ν (Ω) and then s(t) ∈ C ν (Ω).

Consequently Au(t) ∈ C ν (Ω) for t ∈ (0, T ). Remark now that -A e (A i + A e ) -1 A i [u(t)] = Au(t) (with the previous notations). Lemma 21 ensures that A -1 e Au(t) ∈ C 2+ν (Ω), and therefore

(A i + A e )A -1 e Au(t) ∈ C ν (Ω) and at last [u(t)] = -A -1 i (A i + A e )A -1 e Au(t) ∈ C 2+ν (Ω). This implies that x → u(t, x) ∈ C 2 (Ω) for t ∈ (0, T ), since functions in C 2+ν (Ω) belongs to C 2+ν ( Ω).

Global Solution based on a Variational Formulation

The existence of weak solutions for ( 7)-( 12) is established by a Faedo-Galerkin technique: construction of an approximate solution, a priori estimates and compactness result [START_REF] Lions | Résolution de quelques problèmes aux limites non-linéaires[END_REF].

Specific Assumptions and notations

The existence of a weak solution holds under minimal regularity assumptions: Ω has a Lipschitz boundary ∂Ω, σ i,e have L ∞ (Ω) coefficients, and s i,e : [0, +∞) → V ′ are defined and such that s i (t) + s e (t), 1 = 0 for a.e. t > 0, in order to use the bilinear form a and the linear source term s : t ∈ [0, +∞) → s(t) ∈ V ′ as in definition 5.

For sake of simplicity, we assume that m = 1, meaning that w(t, x) ∈ R. We still use the notations

V = H 1 (Ω), H = L 2 (Ω) and U = V /R.
We want to write [START_REF] Hooks | Cardiac microstructure: Implications for electrical propagation and defibrillation in the heart[END_REF] in V ′ and ( 9) in H ′ ≡ H. Hence we need assumptions on f, g : R 2 → R to prove that (u, w) ∈ V × H → (f (u, w), g(u, w)) ∈ V ′ × H ′ is well-defined. Therefore, we shall assume that f and g are both polynomial functions (see lemma 25 below). Additionally, some energy-like estimates are needed to construct the a priori bounds. At last, the technical assumption (H2) below is used to pass to the limit in the variational formulation. We suppose that there exists p ≥ 2 such that the Sobolev embedding [START_REF] Brezis | Analyse fonctionnelle, Théorie and applications[END_REF]; (H2) the functions f and g are affine with respect to w:

V = H 1 (Ω) ⊂ L p (Ω) holds: p ≥ 2 if d = 2; or 2 ≤ p ≤ 6 if d = 3
f (u, w) = f 1 (u) + f 2 (u)w, g(u, w) = g 1 (u) + g 2 w, (48)
where f 1 : R → R, f 2 : R → R, g 1 : R → R are continuous functions and g 2 ∈ R; (H3) there exists constants c i ≥ 0 (i = 1 . . . 6) such that for any u ∈ R,

|f 1 (u)| ≤ c 1 + c 2 |u| p-1 , ( 49
)
|f 2 (u)| ≤ c 3 + c 4 |u| p/2-1 , ( 50
)
|g 1 (u)| ≤ c 5 + c 6 |u| p/2 ; ( 51
) (H4) there exists constants a, λ > 0, b, c ≥ 0 such that for any (u, w) ∈ R 2 , λuf (u, w) + wg(u, w) ≥ a|u| p -b λ|u| 2 + |w| 2 -c. ( 52
)
Remark 24 Three examples of models in electrocardiology satisfying these assumptions will be given in §6.

Using hypothesis (H1), we have the framework

V ⊂ L p (Ω) ⊂ H ≡ H ′ ⊂ L p ′ (Ω) ⊂ V ′ , (53)
meaning in particular that an element u ∈ H

′ or u ∈ (L p (Ω)) ′ is identified to an element of u ∈ H or u ∈ L p ′ (Ω) by u, v = Ω uv.
At last, we use the classical spaces L q (0, T ; X) (1 ≤ q ≤ ∞) of measurable vector valued functions f : t ∈ (0, T ) → f (t) ∈ X where X is a separable Banach space (X alternatively is U, U ′ , V, V ′ or H here). The derivative ∂ t f (or df dt ) of this function is taken in the space of vector valued distributions from (0, T ) onto X [START_REF] Lions | Résolution de quelques problèmes aux limites non-linéaires[END_REF]. A distribution f and a function f ∈ L q (0, T ; X) are identified if

f, φ = T 0 f (t)φ(t)dt (in X) ∀φ ∈ D(0, T ),
where D(0, T ) is the space of real functions C ∞ on R with compact support in (0, T ). We recall that if f ∈ L 1 (0, T ; X) and ∂ t f ∈ L 1 (0, T ; X), then f is equal a.e. to a function in C 0 ([0, T ], X).

We have the Lemma 25 Under hypotheses (H2) and (H3), the mappings

(u, v) ∈ L p (Ω) × H → f (u, v) ∈ L p ′ (Ω) and (u, v) ∈ L p (Ω) × H → g(u, v) ∈ H are well defined.
Specifically, for any (u, w) ∈ L p (Ω) × H, we have

f (u, w) L p ′ (Ω) ≤ A 1 |Ω| 1/p ′ + A 2 u p/p ′ L p (Ω) + A 3 w 2/p ′ H , (54)
g(u, w) H ≤ B 1 |Ω| 1/2 + B 2 u p/2 L p (Ω) + B 3 w H , (55)
where the A i ≥ 0 (i = 1 . . . If p = 2 it is proved by the inequality of Young that

|w||u| p/2-1 ≤ |w| β β + |u| (p/2-1)β ′ β ′ ,
where β = 2/p ′ > 1 and

1 β + 1 β ′ = 1. Since p 2 -1 β ′ = p 2 -1 2 p-1 p-2 = p -1, we have |f (u, w)| ≤ c 1 + c 2 + c 4 β ′ |u| p-1 + c 3 |w| + c 4 β |w| β .
But β > 1 and then we also have |w| ≤ |w| β β + 1 β ′ , so that it can be found positive constants A 1 , A 2 and A 3 such that

|f (u, w)| ≤ A 1 + A 2 |u| p-1 + A 3 |w| β . If p = 2, this inequality is still valid, with A 1 = c 1 , A 2 = c 2 , A 3 = c 3 + c 4 . Now for (u, v) ∈ L p (Ω) × H, we can write f (u, w) L p ′ (Ω) ≤ A 1 + A 2 |u| p-1 + A 3 |w| β L p ′ (Ω) ≤ A 1 L p ′ (Ω) + A 2 |u| p-1 L p ′ (Ω) + A 3 |w| β L p ′ (Ω) = A 1 |Ω| 1/p ′ + A 2 u p/p ′ L p (Ω) + A 3 w 2/p ′ H , (56)
because (p -1)p ′ = p, βp ′ = 2, and similarly,

g(u, w) H ≤ B 1 + B 2 |u| p/2 + B 3 |w| H ≤ B 1 H + B 2 |u| p/2 H + B 3 |w| H = B 1 |Ω| 1/2 + B 2 u p/2 L p (Ω) + B 3 w H . (57)

Existence for the Initial Value Problem

Under the minimal regularity assumptions on the data Ω, σ i,e and s i,e given at the beginning of §5, we are able to write the Definition 26 (Weak solutions) Consider τ > 0 and the three functions

u : t ∈ [0, τ ) → u(t) ∈ H, u e : t ∈ [0, τ ) → u e (t) ∈ H, w : t ∈ [0, τ ) → w(t) ∈ H.
Given (u 0 , w 0 ) ∈ H, we say that (u, u e , w) is a weak solution to (7)-( 12) iff, for any T ∈ (0, τ ),

(1) u : [0, T ] → H and w : [0, T ] → H are continuous, and u(0) = u 0 , w(0) = w 0 in H (that is eq. ( 12)); (2) for a.e. t ∈ (0, τ ), we have u(t) ∈ V , u e (t) ∈ V /R, and u ∈ L 2 (0, T ;

V ) ∩ L p (Q T ), where Q T = (0, T ) × Ω;
and (u, u e , w) verify in D ′ (0, T ):

d dt (u(t), v) + Ω σ i ∇(u(t) + u e (t)) • ∇v + Ω f (u(t), w(t))v = s i (t), v , d dt (w(t), v) + Ω g(u(t), w(t))v = 0,
respectively for all v ∈ V and for all v ∈ H, and

Ω σ i ∇u(t) • ∇v e + Ω (σ i + σ e )∇u e (t) • ∇v e = s i (t) + s e (t), v e , ∀v e ∈ V /R. (58)
Remark 27 The weak derivatives of u :

t ∈ [0, T ] → H and w : t ∈ [0, T ] → H identify to functions ∂ t u ∈ L 2 (0, T ; V ′ ) + L p ′ (Q T ) and ∂ t w ∈ L 2 (0, T ; V ′ ) + L p ′ (Q T ).
Indeed the following equalities are true in D ′ (0, T ):

∂ t u, v = d dt (u(t), v) ∀v ∈ V = V ∩ L p (Ω), ∂ t w, v = d dt (w(t), v) ∀v ∈ H.
Naturally, we have the two following lemma:

Lemma 28 The functions (u, u e , w) are a weak solution to (7)-( 12) iff conditions ( 1)-( 2) of definition 26 hold and (u, w) verify in D ′ (0, T ):

d dt (u(t), v) + a(u(t), v) + Ω f (u(t), w(t))v = s(t), v ∀v ∈ V, d dt (w(t), v) + Ω g(u(t), w(t))v = 0 ∀v ∈ H,
where a(•, •) and s ∈ V ′ are given in definition 5. The function u e is then recovered from (58).

Lemma 29 Any strong solution (u, u e , w) on [0, τ ) is a weak solution on [0, τ ). Conversely, if ∂Ω is C 1 regular, any weak solution (u, u e , w) on [0, τ ) such that u(t) ∈ H 2 (Ω) for a.e. t ∈ [0, τ ) is a strong solution.
Theorem 30 (Global existence of a weak solution) Let Ω, σ i,e have the minimal regularity specified in §2. Suppose that hypotheses (H1) to (H4) on f, g hold for some p ≥ 2. Let be given u 0 , w 0 ∈ H and s i , s e ∈ L 2 (R + ; V ′ ) such that s i (t) + s e (t), 1 for a.e. t > 0. Then the system (7)-( 12) has a weak solution (u, u e , w) in the sense of definition 26 with τ = +∞.

Proof. Using lemma 28, it is given in the next subsections, in three parts:

• construction of an approximate solution using the Faedo-Galerkin technique; • a priori estimates on the approximate solution;

• compactness results, and convergence of the approximate solution towards a weak solution.

Construction of an approximate solution

We use the special orthonormal Hilbert basis (in H) (ψ i ) i∈N of eigenvectors defined in theorem 6. For m ≥ 1, we note V m = span(ψ 0 , . . . , ψ m ) ⊂ V . We are looking for a couple of functions t → (u m (t), w m (t)) with

u m (t) = m i=0 u im (t)ψ i ∈ V m , w m (t) = m i=0 w im (t)ψ i ∈ V m
where (u im (t), w im (t)) i=0...m are real valued functions solutions of

d dt u im (t) + λ i u im (t) + Ω f (u m (t), w m (t))ψ i = s(t), ψ i , (59)
d dt w im (t) + Ω g(u m (t), w m (t))ψ i = 0 (60)
for i = 0 . . . m, and with initial data

u m (0) = u m0 , w m (0) = w m0 . (61)
Since u 0 and w 0 are in H, we can take u m0 and w m0 to be the H orthogonal projections of u 0 and w 0 on V m :

u m0 -u 0 H → 0, w m0 -w 0 H → 0 as m → ∞. (62)
Equations ( 59) and ( 60) make sense because

u m (t) ∈ V ⊂ L p (Ω), w m (t) ∈ H so that f (u m (t), w m (t)) ∈ L p ′ (Ω) ⊂ V ′ and g(u m (t), w m (t)) ∈ H (from lemma 25
) and ψ i ∈ V ⊂ L p (Ω) and ψ i ∈ H. Moreover, it can easily be seen that the last three terms in Eq. ( 59) and the last term in Eq. ( 60) are continuous functions of u im and w im .

The initial value problem composed of the 2m + 2 differential equations ( 59)-(60) with initial data (61) has a maximal solution defined for t ∈ [0, t m ) with u im and w im in C 1 (theorem of Cauchy-Peano from [43, p.59]). If (u m , w m ) is not a global solution (i.e. t m < +∞) then it is unbounded in [0, t m ). It will be shown in the next section, using a priori estimates, that (u m , w m ) remains bounded for all time, namely t m = +∞.

A priori estimates

The following lemma establishes uniform bounds for any T ∈ (0, t m ), first on the sequences u m and w m in L ∞ (0, T ; H), then on the sequences u m , u ′ m respectively in L p (Q T ) ∩ L 2 (0, T ; V ) and its dual L p ′ (Q T ) + L 2 (0, T ; V ′ ), and finally on the sequences w m , w ′ m in L 2 (Q T ) (identified to its dual). We use the norm

• L p (Q T )∩L 2 (0,T ;V ) = max( • L p (Q T ) , • L 2 (0,T ;V ) ) and the dual norm u L p ′ (Q T )+L 2 (0,T ;V ′ ) = inf u=u 1 +u 2 u 1 L p ′ (Q T ) + u 2 L 2 (0,T ;V ′ ) .
Lemma 31 (A priori estimates) The maximal solution of the Cauchy problem (59)-( 61) is defined for any t > 0; and for any T > 0, there exists positive

constants C 1 , C 2 , C 3 , C 4 such that λ u m (t) 2 H + w m (t) 2 H ≤ C 1 , ∀t ∈ [0, T ], (63) u m L p (Q T )∩L 2 (0,T ;V ) ≤ C 2 , (64) u ′ m L p ′ (Q T )+L 2 (0,T ;V ′ ) ≤ C 3 , (65)
w ′ m L 2 (Q T ) ≤ C 4 , (66)
where

u ′ m (t) = m i=0 u ′ im (t)ψ i , w ′ m (t) = m i=0 w ′ im (t)ψ i are the derivative of u m : [0, +∞) → V and w m : [0, +∞) → H.
The estimate (63) is the bound in L ∞ (0, T ; H) for u m and w m ; and the bound for w m in L 2 (Q T ) is easily derived from it.

Proof. Multiplying (59) by λu im (λ > 0 defined in hypothesis (H4)), multiplying (60) by w im , and summing over i = 1 . . . m yields, for any t ∈ [0, t m ),

1 2 λ d dt u m 2 H + 1 2 d dt w m 2 H + λa(u m , u m ) + Ω (λf (u m , w m )u m + g(u m , w m )w m ) = λ s, u m .
Using the properties of a(•, •) from theorem 6 (coercivity and continuity) and hypothesis (H4), we have for any t ∈ [0, t m ) and for any ξ > 0, 1 2

d dt λ u m (t) 2 H + w m (t) 2 H + αλ u m (t) 2 V + a Ω |u m (t)| p ≤ (b + α) λ u m (t) 2 H + w m (t) 2 H + c|Ω| + s(t) V ′ u m (t) V ≤ (b + α) λ u m (t) 2 H + w m (t) 2 H + c|Ω| + 1 2ξ s(t) 2 V ′ + ξ 2 u m (t) 2 V .
And then, choosing ξ = αλ, we have 1 2

d dt λ u m (t) 2 H + w m (t) 2 H + αλ 2 u m (t) 2 V + a Ω |u m (t)| p ≤ b λ u m (t) 2 H + w m (t) 2 H + c|Ω| + 1 2αλ s(t) 2 V ′ , (67)
with b = b + α.

We know that u m (0) H ≤ u 0 H , w m (0) H ≤ w 0 H , Ω is bounded and S t := t 0 s(τ ) 2 V ′ dτ < +∞. From the Gronwall inequality there exists a constant C 1 > 0 that depends only on σ i,e , f , g, u 0 , w 0 , Ω, s i,e and t m , such that 0

≤ t < t m ⇒ λ u m (t) 2 H + w m (t) 2 H ≤ C 1 .
As an immediate consequence the solution u m , w m cannot explode in finite time. It is defined on [0, +∞) (it is a global solution, t m = +∞). Now, for any fixed T > 0, we have found a constant C 1 > 0 such that (63) is valid. Coming back with C 1 into (67) we immediately have the estimate (64) of lemma 31 with

C 2 = max 2C T αλ 1/2 , C T a 1/p , where C T = 1 2 λ u 0 2 H + w 0 2 H + bT C 1 + cT |Ω| + 1 2αλ S T .
The remaining estimates on u ′ m , w ′ m is the most difficult. Consider the projection P m : V ′ → V ′ defined for u ∈ V ′ by

P m u = m i=1 u, ψ i ψ i .
It is equivalently defined as the unique element of V m such that u, v = P m u, v for all v ∈ V m . For any v ∈ V and any t > 0, remark that

d dt (u m (t), v) = (u ′ m (t), v) = u ′ m (t), v , Ω f (u m (t), w m (t))v = f (u m (t), w m (t)), v , because u ′ m (t) ∈ V m ⊂ V ′ and f (u m (t), w m (t) ∈ L p ′ (Ω) while v ∈ V ⊂ L p (Ω). And then equation (59) reads ∀v ∈ V m , ∀t > 0, u ′ m (t), v = -Au m (t) + f (u m (t), w m ), v + s(t), v , so that ∀t > 0, u ′ m (t) = -P m (Au m (t) + f (u m (t), w m (t)) + s(t)) , ( 68
)
where A is the weak operator associated to the bilinear form a(•, •) on V × V , as defined in lemma 9.

For T > 0 fixed, we have from (64) and the continuity of A,

Au m L 2 (0,T ;V ′ ) ≤ M T 0 u m (t) 2 V dt 1/2 ≤ MC 2
and from (63), (64) and lemma 25,

f (u m , w m ) L p ′ (Q T ) ≤ A 1 |Ω| 1/p ′ + A 2 u m (t) p/p ′ L p (Ω) + A 3 w m (t) 2/p ′ H L p ′ (0,T ) ≤ A 1 (|Ω|T ) 1/p ′ + A 2 u m p/p ′ L p (Q T ) + A 3 w m 2/p ′ L 2 (Q T ) ≤ A 1 (|Ω|T ) 1/p ′ + A 2 C p/p ′ 2 + A 3 (C 1 T ) 1/p ′ .
It remains to bound the projection operator P m . First, remark that the restriction of P m to V can be viewed as an operator from V onto V (since

P m (V ′ ) ⊂ V m ⊂ V ), given by ∀u ∈ V, P m u = m i=1 (u, ψ i )ψ i .
For u ∈ H, P m u is the orthogonal projection of u on V m , and P m u H ≤ u H . The transpose P T m of P m|V identifies with P m : V ′ → V ′ (simple computation), and then we have

P m L(V ′ ,V ′ ) = P m L(V,V) . For u ∈ V we can compute a(P m u, P m u) = +∞ i=0 λ i (P m u, ψ i )(P m u, ψ i ) = m i=0 λ i (u, ψ i )(u, ψ i ) ≤ +∞ i=0 λ i |(u, ψ i )| 2 = a(u, u).
As a consequence, for all u ∈ V ,

α P m u 2 V ≤ a(P m u, P m u) + α P m u 2 H ≤ M u 2 V + α u 2 H ≤ (M + α) u 2 V .
It shows that P m is uniformly bounded in V ′ : P m L(V ′ ,V ′ ) ≤ 1 + M α , and we have

P m (Au m ) L 2 (0,T ;V ′ ) ≤ 1 + M α MC 2 , P m (f (u m , w m )) L p ′ (Q T ) ≤ 1 + M α A 1 (|Ω|T ) 1/p ′ + A 2 C p/p ′ 2 + A 3 (C 1 T ) 1/p ′ , P m s L 2 (0,T ;V ′ ) ≤ 1 + M α s L 2 (0,T ;V ′ )
The bound (65) is a consequence of these estimates and of (68).

In a similar way, equation (60) reads

∀v ∈ V m ⊂ H, ∀t > 0, w ′ m (t), v = -g(u m (t), w m (t)), v , so that ∀t > 0, w ′ m (t) = -P m (g(u m (t), w m (t))) , (69)
where the operator P m can be restricted to the orthogonal projection P m|H , in particular, P m L(H,H) ≤ 1.

For T > 0 fixed, from (63), (64) and lemma 25, we have (66):

w ′ m L 2 (Q T ) ≤ g(u m , w m ) L 2 (Q T ) ≤ B 1 |Ω| 1/2 + B 2 u m (t) p/2 L p (Ω) + B 3 w m (t) H L 2 (0,T ) ≤ B 1 (|Ω|T ) 1/2 + B 2 u m p/2 L p (Q T ) + B 3 w m L 2 (Q T ) ≤ B 1 (|Ω|T ) 1/2 + B 2 (C 2 ) p/2 + B 3 (C 1 T ) 1/2 := C 4 .

Convergence towards a solution

It is easy to see that

L p ′ (Q T ) + L 2 (0, T ; V ′ ) ⊂ L p ′ (0, T ; V ′ ) since p ′ ≤ 2 and L p ′ (Ω) ⊂ V ′ (Sobolev inequality).
Hence the sequence (u ′ m ) remains in a bounded set of L p ′ (0, T ; V ′ ) while (u m ) remains in a bounded set of L 2 (0, T ; V ). It follows from a classical compactness result, see for instance [39, th. 5.1, p.58], that the sequence (u m ) has a subsequence that converges in L 2 (Q T ).

As a consequence, we can construct subsequences of u m and w m , still denoted by u m and w m , such that

• u m → u weak in L p (Q T ) ∩ L 2 (0, T, V ) and u ′ m → ũ weak in L p ′ (Q T ) + L 2 (0, T, V ′ ), • w m → w weak in L 2 (Q T ), and w ′ m → w weak in L 2 (Q T ),
and from the compactness result,

• u m → u strong in L 2 (Q T ), and then almost everywhere in Q T ,

where u ∈ L p (Q T ) ∩ L 2 (0, T, V ), w ∈ L 2 (Q T ), and ũ ∈ L p ′ (Q T ) + L 2 (0, T, V ′ ), w ∈ L 2 (Q T ).
For i ≥ 1 fixed and φ ∈ D(0, T ), we naturally have

- T 0 Ω u ′ m ψ i φ = T 0 Ω u m ψ i φ ′ → T 0 Ω uψ i φ ′ , - T 0 Ω w ′ m ψ i φ = T 0 Ω w m ψ i φ ′ → T 0 Ω wψ i φ ′ , because ψ i φ ′ ∈ L 2 (Q T ) ∩ L p (Q T ) ∩ L 2 (0, T ; V ).
As a consequence, we have in the space D ′ (0, T ) of distribution on (0, T ),

d dt (u(t), ψ i ) = ũ(t), ψ i , d dt (w(t), ψ i ) = ( w(t), ψ i ) . (70)
Since a(•, •) is bilinear and continuous on V ×V and ψ i φ ∈ L p (Q T )∩L 2 (0, T ; V ) for any φ ∈ D(0, T ), we have

∀φ ∈ D(0, T ), T 0 a(u m (t), φ(t)ψ i )dt → T 0 a(u(t), φ(t)ψ i )dt.
Concerning the nonlinear terms, we use hypothesis (H2) to write

f (u m , w m ) = f 1 (u m ) + f 2 (u m )w m = f 1 (u m ) + (f 2 (u m ) -f 2 (u)) w m + f 2 (u)w m , g(u m , w m ) = g 1 (u m ) + g 2 w m . Now, we have u m → u a.e. in Q T and f 1 is continuous, so that f 1 (u m ) → f 1 (u) a.e.in Q T ; and f 1 (u m ) is uniformly bounded in L p ′ (Q T ), f 1 (u m ) L p ′ (Q T ) ≤ c 1 + c 2 |u m | p-1 L p ′ (Q T ) ≤ c 1 (|Ω|T ) 1/p ′ + c 2 u m p/p ′ L p (Q T ) .
It follows from a classical result, see [39, lemma 1.3, p.12], that f

1 (u m ) → f 1 (u) weak in L p ′ (Q T ): ∀φ ∈ D(0, T ), T 0 (f 1 (u m (t)), φ(t)ψ i )dt → T 0 (f 1 (u(t)), φ(t)ψ i )dt.
Similarly, g 1 is continuous and then g 1 (u m ) → g 1 (u) a.e.in Q T ; and g

1 (u m ) is bounded in L 2 (Q T ), g 1 (u m ) L 2 (Q T ) ≤ c 5 + c 6 |u m | p/2 L 2 (Q T ) ≤ c 5 (|Ω|T ) 1/2 + c 6 u m p/2 L p (Q T ) ,
and then g

1 (u m ) → g 1 (u) weak in L 2 (Q T ), ∀φ ∈ D(0, T ), T 0 (g 1 (u m (t)), φ(t)ψ i )dt → T 0 (g 1 (u(t)), φ(t)ψ i )dt. Since w m → w weak in L 2 (Q T ) we naturally have ∀φ ∈ D(0, T ), T 0 (g 2 w m (t), φ(t)ψ i )dt → T 0 (g 2 w(t), φ(t)ψ i )dt. As f 2 (u)φ(t)ψ i ∈ L 2 (Q T ) from hypothesis (H3), the weak convergence of w m in L 2 (Q T ) also implies that ∀φ ∈ D(0, T ), T 0 (f 2 (u)w m (t), φ(t)ψ i )dt → T 0 (f 2 (u)w(t), φ(t)ψ i )dt.
But we also have the following identity in D ′ (0, T ):

∂ t u(t), u(t) = 1 2 d dt u(t) 2 H ,
and then 1 2

d dt u(t) 2 H = -a(u(t), u(t)) -f (u(t), w(t)), u(t) + s(t), u(t) , so that t → u(t) 2
H is H 1 (0, T ), and then it is continuous from [0, T ] to R. As a consequence, the function u : t ∈ [0, T ] → u(t) ∈ H is continuous. Since u m (0) → u 0 and w m (0) → w 0 in H, we easily prove that u(0) = u 0 and w(0) = w 0 .

Uniqueness

Assume that (u 1 , u e1 , w 1 ) and (u 2 , u e2 , w 2 ) are two weak solutions of ( 7)-( 12) with the same initial data u 1 (0) = u 2 (0) = u 0 and w

1 (0) = w 2 (0) = w 0 . For any u ∈ L 2 (0, T ; V ) ∩ L p (Q T ) and w ∈ L 2 (Q T ) with ∂ t u ∈ L 2 (0, T ; V ′ ) + L p ′ (Q T ) and ∂ t w ∈ L 2 (Q T ), we have in D ′ (0, T ) that ∂ t u(t), u(t) = 1 2 d dt u(t) 2 H , ∂ t w(t), w(t) = 1 2 d dt w(t) 2 H .
As a consequence, we easily prove that 1 2

d dt u 1 -u 2 2 H + a(u 1 -u 2 , u 1 -u 2 ) + Ω (f (u 1 , w 1 ) -f (u 2 , w 2 )) (u 1 -u 2 ) = 0, and 1 2 d dt w 1 -w 2 2 H + Ω (g(u 1 , w 1 ) -g(u 2 , w 2 )) (w 1 -w 2 ) = 0.
With a linear combination of these two equations, we will be able to conclude using a Gronwall inequality if we can bound below Φ(u 1 , w 1 , u 2 , w 2 ) = Now, let Q(z) = 1 2 (∇F (z) T +∇F (z)) be the symmetric part of ∇F for z ∈ R 2 , and denote by λ 1 (z) ≤ λ 2 (z) its eigenvalues. We can complete the proof under the hypothesis that

∃C ∈ R, ∀z ∈ R 2 , λ 2 (z) ≥ λ 1 (z) ≥ C. (71)
Indeed, in that case, we have for any

z 1 , z 2 ∈ R 2 , Φ(z 1 , z 2 ) = Ω 1 0 (z 1 -z 2 ) T [∇F (z θ )](z 1 -z 2 ) dθdx ≥ C Ω 1 0 |z 1 -z 2 | 2 dθdx ≥ C min(1, µ -1 ) µ u 1 -u 2 2 H + w 1 -w 2 2 H .
As a consequence, taking Y (t) = (µ u 1 (t) -u 2 (t) 2 H + w 1 (t) -w 2 (t) 2 H ), under assumption (71) on the data f and g, we prove that Using the lemma of Gronwall, we have proved the following result:

Theorem 32 If the condition (71) is satisfied, then the solution obtained in Theorem 30 is unique.

Remark that Eq. (72) also provide a stability estimate with respect to the initial condition.

We will apply that theorem to the first example presented below, though it is not clear how to obtain uniqueness for the last two examples.

Examples

FitzHugh-Nagumo

The FitzHugh-Nagumo model reads as f (u, w) = u(u -a)(u -1) + w, g(u, w) = -ǫ(ku -w), with 0 < a < 1, k, ǫ > 0. The functions f and g are obviously of the form (48) with f 1 , f 2 , g 1 continuous and g 2 = ǫ. Using Young's inequality, we have Here, we will prove (52) because it allows negative |u| 2 and |w| 2 bounds below, so that the terms in u 3 and u 2 w can be manipulated to enter the main bound λku 4 . For instance with λ = ǫk, we write

|u| 2 ≤ 2|u| 3 3 + 1 3 , |u| ≤ |u| 3 3 + 2 3 , |u| ≤ |u| 2 2 + 1 2 , (73)
|

3 )

 3 and B i ≥ 0 (i = 1 . . . 3) are numerical constants that depend only on the c i (i = 1 . . . 6) and on p. Proof. For (u, w) ∈ R 2 , we have from (H2) and (H3), |f (u, w)| ≤ c 1 + c 2 |u| p-1 + c 3 |w| + c 4 |w||u| p/2-1 , |g(u, w)| ≤ B 1 + B 2 |u| p/2 + B 3 |w|, with exactly B 1 = c 5 , B 2 = c 6 and B 3 = |g 2 |.

Ω µ (f (u 1 , w 1 )F 2 F (z 1 ) 1 0[

 11211 -f (u 2 , w 2 )) (u 1 -u 2 ) + (g(u 1 , w 1 ) -g(u 2 , w 2 )) (w 1 -w 2 )dx for some µ > 0. Consider the function F : R 2 → R 2 defined by (u, w) ∈ R 2 ,and denote by z ∈ R 2 the vector z = (u, w) T ∈ R 2 . Then we haveΦ(u 1 , w 1 , u 2 , w 2 ) = Φ(z 1 , z 2 ) = Ω (F (z 1 ) -F (z 2 )) • (z 1 -z 2 )dx,where • denotes the inner product in R 2 . Here F is continuously differentiable, so that Taylor expansion with an integral remainder implies that ∀z 1 , z 2 ∈ R -F (z 2 ) = ∇F (z θ )](z 1 -z 2 ) dθ where z θ = θz 1 + (1 -θ)z 2 and ∇F =

1 2 Y

 12 ′ (t) ≤ -C min(1, µ -1 )Y (t),(72)for any t ∈ [0, T ].

2 .|u| 4 ,

 24 and then (H3) holds with p = 4 (and c 4 = 0):|f 1 (u)| = |u(u -a)(u -1)| ≤ 2 a) + 1 |u| 3 , |f 2 (u)| = 1, |g 1 (u)| = ǫk|u| ≤ Consider the function E(u, w) = ǫkuf (u, w)+wg(u, w) defined in R 2 .We haveE(u, w) = ǫku 4 -ǫk(1 + a)u 3 + ǫkau 2 + ǫw 2 ≥ ǫk |u| 4 -(1 + a)|u| 3 .With Young's inequality, we can find a constant γ > 0 such that(1 + a)|u| 3 ≤ |u| 4 2 + γ.Consequently,E(u, w) + ǫkγ ≥ ǫk 2 which is exactly (H4) with λ = kǫ, a = kǫ/2, b = 0 and c = kǫγ. As regards the uniqueness of the solution, we verify the condition (71) to apply Theorem 32. One easily calculates ∇F (Taking µ = ǫk, we get rid of the antisymmetric part in the quadratic form and easily bound below the eigenvalues by C = ǫ min(k(a -(1 + a) 2 /3), 1).

6. 2

 2 Aliev-PanfilovThe Aliev-Panfilov model[START_REF] Panfilov | A simple two-variable model of cardiac excitation[END_REF] isf (u, w) = ku(u -a)(u -1) + uw, g(u, w) = ǫ (ku(u -1 -a) + w) ,with 0 < a < 1, k, ǫ > 0. The functions f and g are obviously of the form (48) with f 1 , f 2 , g 1 continuous and g 2 = ǫ. Using the inequalities (73), we get (H3) with p = 4 (and c 4 = 1, c 3 = 0):|f 1 (u)| = k|u(u -a)(ua) + 1 k|u| 3 , |f 2 (u)| = |u|, |g 1 (u)| = ǫk|u(u -1 -a)| ≤ 1 2 ǫk(1 + a) + 1 2 (1 + a) + 1 ǫk|u| 2 .Now, we compute the function E(u, w) = λuf (u, w) + wg(u, w). It is E(u, w) = λku 4 -λk(1+a)u 3 +λkau 2 +(λ+ǫk)u 2 w -ǫk(1+a)uw +ǫw 2 . (74)

  ×B m ; and the source term S : t ∈ [0, +∞) → (s(t), 0) ∈ Z where A and s(t) are given in definition 12. Since the operators A and w ∈ B m → 0 ∈ B m are self-adjoint, they are also sectorial. Thus A is sectorial as the Cartesian product of two sectorial operators. See[31, p.18-19] for the definition and properties of sectorial operators.

	Lemma 14 The unbounded operator A : D(A) ⊂ H → H is a sectorial
	operator.
	Proof.

Lemma 15 If s i,e : [0, +∞) → H are locally ν-Hölder continuous functions with s i (t)+s e (t) ∈ H/R for all t ≥ 0, then S : [0, +∞) → Z is locally ν-Hölder continuous. Proof. Consider [t 1 , t 2 ] ⊂ [0, +∞), and the constant C > 0 such that (45) holds. If δ 1 , δ 2 ∈ [t 1 , t 2 ], we have

  + ǫ|w| 2 + ǫk 2 a|u| 2 . = ǫkβ 2 , second authors, respectively to the INRIA-Rocquencourt and the University of Ottawa, and the Centre de recherches mathématiques (CRM) of the Université de Montréal for providing a post-doctoral fellowships to the third author. This work was also supported by a NSERC Discovery Grant.

	(1 + a)u 3 | ≤ |u 2 w| ≤	3 4 1 2	α|u| 3 4/3 + (β|u| 2 ) 2 + 1 2	1 4 |w| 1 + a α 2 β ,	4	,	(75) (76)
	|uw| ≤	1 2	|u| 2 +	1 2	|w| 2 ,	(77)
	for any α > 0 and β > 0, and then
	E(u, w) ≥ ǫk 2 -ǫk 2 3 4 -1 4 ǫk 2 1 + a α 4 -ǫk 1 β 2 |w| 2 -ǫk α 4/3 -ǫkβ 2 |u| 4 1 + a 2 |u| 2 -ǫk |w| 2 Now, we can simply take α and β such that 1 + a 2
	3 4	α 4/3 =	1 2	, and	1 4	ǫk 2

∃C > 0, ∀x, y ∈ Ω, |w(y)w(x)| ≤ C|x -y| ν .

Acknowledgments

The authors would like to thank the Laboratoire Jean Leray, Université de Nantes, for supporting a one-month stay of the first author in Nantes, the INRIA associated team ACE who covered travel expenses of the first and

⋆ This work was supported by a NSERC Discovery Grant and the INRIA associated team ACE.

tinuous on R × R m , we have

and in that case, F : z ∈ Z α → F (z) ∈ Z is locally Lipschitz continuous.

Lemma 17 (Case B = C ν (Ω)) For B = C ν (Ω), f, g C 2 functions on R × R m , and α < 1, we have

and in that case, F : z ∈ Z α → Z is locally Lipschitz continuous.

Proof. From lemma 14 the operatorA is a sectorial operator, hence the regularity embeddings from [31, p.39] apply.

A locally Lipschitzian function f : R → R induces a locally Lipschitzian function f : L ∞ (Ω) → L ∞ (Ω), so that F can be extended to a locally Lipschitz continuous function

. The mapping F can be extended to a locally Lipschitzian function F : Z α → Z.

Existence of and uniqueness of local in time solution

We are ready to define the notion of strong solution local in time:

Definition 18 (Strong solution) Consider τ > 0 and the functions z : t ∈ [0, τ ) → z(t) = (u(t), w(t)) ∈ Z and u e : t ∈ [0, τ ) → u e (t) ∈ H. Given (u 0 , w 0 ) ∈ Z, we say that (u, u e , w) is a strong solution to (7)-( 12) iff, and(u, u e , w) verify (7)-( 9) for all t ∈ (0, τ ) and for a.e. x ∈ Ω, and the boundary conditions (10)-( 11) for all t ∈ (0, τ ) and for a.e. x ∈ ∂Ω.

One easily derives the following characterization

Lemma 19

The functions z = (u, w) and u e are a strong solution to (7)- [START_REF]Computational Biology of the Heart[END_REF] iff conditions (1)-(3) of definition 18 and condition (4') below are satisfied:

The remaining term in f is such that

Remark that

The duality product on the right hand side makes sense because

and we can conclude with the same classical result that (f

Gathering all these results, the functions u and w verify, for any i ≥ 1,

in the space of distributions D ′ (0, T ), for any i ≥ 0. Since (ψ i ) i≥0 is dense in V , this is exactly the desired result (lemma 28).

Continuity

We have

; and we deduce from (70) that u and w have their weak derivative ∂ t u and ∂ t w respectively in L 2 (0, T ;

It is deduced from a classical result, see for instance [39,

MacCulloch

The model introduced by McCulloch [START_REF] Roger | A collocation-Galerkin finite element model of cardiac action potential propagation[END_REF] is

with 0 < a < 1, and b, k, ǫ > 0. The functions f and g are obviously of the form (48) with f 1 , f 2 , g 1 continuous and g 2 = ǫ. Using the inequalities (73), we get (H3) with p = 4 (and c 4 = 1, c 3 = 0):