Mostafa Bendahmane

Fatima Mroue

AND Mazen Saad

Raafat Talhouk

Unfolding homogenization method applied to physiological and phenomenological bidomain models in electrocardiology

Keywords: Mathematics Subject Classification. 65N55, 65M Bidomain model, reaction-diffusion system, Homogenization theory, unfolding method, convergence

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The heart is the muscular organ that contracts to pump blood throughout the body. Its contraction is initiated by an electrical signal called action potential. At a microscopic level, the cardiac tissue is a complex structure composed of elongated connected cells (cardiomyocytes) that have a cylindrical shape and that are aligned in preferential directions forming fibers. Cardiomyocytes are encapsulated in a dynamic cell membrane (the sarcolemma) that separates the interior of the cell from the surrounding medium and maintains a potential difference (the transmembrane potential) between the two media due to the different concentrations of various ionic species on both sides. The elongated cardiomyocytes are endowed with special end-to-end connections (the gap junctions) that form the long fiber structure of the muscle, as well as with lateral junctions that permit the connection between the intracellular spaces of the elongated fibers. Since those connections have a low resistance, the cardiac tissue can be viewed as a single intracellular connected domain, separated from the extracellular domain by the surface of the cell membrane [29]. Moreover, the sarcolemma consists of a phospholipid bilayer in which are embedded ionic channels that ensure the flow of ionic currents from the extra-to intracellular space or vice versa. As a consequence of this transfer of ionic species between the two-spaces (intra-and extracellular spaces) a current flows across the cell membrane (transmembrane current). The capacitive, diffusive and conductive effects contribute to this current flux across the membrane [29,[START_REF] Clayton | A guide to modelling cardiac electrical activity in anatomically detailed ventricles[END_REF][START_REF] Sundnes | Computing the electrical activity in the heart[END_REF]. From a physical point of view, the cardiac tissue can be viewed as partitioned into two ohmic conducting volumes (intra-and extracellular spaces). The intra-and extracellular domains act as volume conductors and can be described by a quasi-static approximation of elliptic equations in both spaces. These equations are complemented by a dynamical boundary equation at the interface of the two regions. It is worth mentioning in the sequel that the approximation of the ionic current flow is based on Ohm's law and charge conservation and that these equations depend (at the microscopic level) on a small parameter (0 < ε << 1) whose order of magnitude is the ratio of the two macro-and microscopic space scales. In this paper we derive a macroscopic bidomain model of cardiac electrophysiology based on a microscopic bidomain model, using a rigorous homogenization method. Indeed, the microscopic model is unsuitable for numerical computations due to the complexity of the underlying geometry, which highlights the importance of the rigourous derivation of the macroscopic model while taking into account the properties of the physiological and microscopic structure. Classically, homogenization has been done by means of the multiple-scale method which permits to formally obtain the homogenized problem based on a formal asymptotic expansion [START_REF] Cioranescu | An introduction to homogenization, volume 17 of oxford lecture series in mathematics and its applications[END_REF][START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF]. There are now various mathematical methods related to this theory: the oscillating test functions method due to L. Tartar in [START_REF] Tartar | Quelques remarques sur l'homogénéisation[END_REF], the two-scale convergence method introduced by G. Nguetseng in [START_REF] Nguetseng | A general convergence result for a functional related to the theory of homogenization[END_REF], and further developed by G. Allaire in [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] (see also [START_REF] Allaire | Homogenization of the neumann problem with nonisolated holes[END_REF]) and recently the periodic unfolding method introduced by D. Cioranescu, A. Damlamian and G. Griso for the study of classical periodic homogenization in the case of fixed domains and adapted to homogenization in domains with holes in [START_REF] Cioranescu | The periodic unfolding method in perforated domains[END_REF]. The idea of the unfolding operator was used in [START_REF] Bourgeat | Convergence of the homogenization process for a double-porosity model of immiscible two-phase flow[END_REF][START_REF] Arbogast | Derivation of the double porosity model of single phase flow via homogenization theory[END_REF][START_REF] Vogt | A homogenization theorem leading to a Volterra-integrodifferential equation for permeation chromatography[END_REF] under the name of periodic modulation or dilation operator. The name "unfolding operator" was then introduced in [START_REF] Cioranescu | The periodic unfolding method in perforated domains[END_REF] and deeply studied in [START_REF] Cioranescu | The periodic unfolding method in homogenization[END_REF][START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF]. The interest of the unfolding method comes, on one hand, from the fact that it only deals with functions and classical notions of convergence in L p spaces and it does not necessitate the use of a special class of test functions. On the other hand, the unfolding operator maps functions defined on oscillating domains into functions defined on fixed domains. Hence, the proof of homogenization results becomes quite simple.

Regarding the asymptotic behavior of a microscopic-level modeling problem for the bioelectric activity of the heart, there is the work by M. Pennachio, G. Savaré, and P. Franzone that rigourously studies the derivation of the bidomain model in the framework of Γ-convergence theory presented in [START_REF] Pennacchio | Multiscale modeling for the bioelectric activity of the heart[END_REF]. Recently, the two-scale method has been used in [START_REF] Collin | Mathematical analysis and 2-scale convergence of a heterogeneous microscopic bidomain model[END_REF][START_REF] Grandelius | The cardiac bidomain model and homogenization[END_REF] to obtain the homogenized macroscopic model using different ionic models and assumptions on the conductivity matrices. In [START_REF] Collin | Mathematical analysis and 2-scale convergence of a heterogeneous microscopic bidomain model[END_REF], the authors derive a macroscopic bidomain model using simplified ionic models whereas in [START_REF] Grandelius | The cardiac bidomain model and homogenization[END_REF], the authors use the FitzHugh-Nagumo ionic model. In the present work, we treat a generalized class of ionic models including the FitzHugh-Nagumo model along with physiological models involving ionic concentrations that appear as arguments of a logarithmic function and that must be shown to be bounded away from 0. We further note that in [START_REF] Collin | Mathematical analysis and 2-scale convergence of a heterogeneous microscopic bidomain model[END_REF][START_REF] Grandelius | The cardiac bidomain model and homogenization[END_REF], the cardiac domain was assumed to be a cube in R 3 . Regarding the mathematical analysis of the microscopic model, we point out that in [START_REF] Veneroni | Reaction-diffusion systems for the microscopic cellular model of the cardiac electric field[END_REF], the author used Schauder's fixed point theorem and in [START_REF] Colli | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF], the authors used a variational approach to establish the well-posedness of the microscopic problem under different initial and boundary conditions. In the present work, we prove the existence of solution of the microscopic problem by a constructive method based on the Faedo-Galerkin approach without the restrictive assumption, usually found in the literature, on the conductivity matrices to have the same basis of eigenvectors or to be diagonal matrices (see for instance [START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF] where the authors prove the existence of a local in time strong solution of the bidomain equations after introducing the so-called bidomain operator). It is worth to mention that our approach is innovative and cannot be found in the literature in the context of existence of solutions to the microscopic bidomain model. The convergence of solutions of a sequence of microscopic problems to the solution of the macroscopic problem is established in properly chosen function spaces. We use the unfolding method in perforated domains [START_REF] Cioranescu | The periodic unfolding method in homogenization[END_REF][START_REF] Cioranescu | The periodic unfolding method in perforated domains[END_REF], for sequences of functions bounded in L 2 , H 1 or in H 1/2 on a micro-periodic domain. The difficulty of the homogenization problem for the bidomain equations is due, on one hand, to the degenerate structure of the equations, in combination with the highly oscillating underlying geometry. As a consequence, standard parabolic a priori estimates are not immediately available [START_REF] Colli | Degenerate evolution systems modeling the cardiac electric field at micro-and macroscopic level[END_REF]. On the other hand, the (nonlinear) dynamics of the cellular model take place on the cell membrane which is a wildly oscillating surface. Hence, an ambiguity arises in defining a proper notion of "strong convergence" of functions in this context. However, some kind of strong convergence is required to pass to the limit in the nonlinear equations. For this reason, we also use the boundary unfolding operator along with a Kolmogorov-Riesz compactness argument [START_REF] Amann | Compact embeddings of vector valued sobolev and besov spaces[END_REF][START_REF] Gahn | A characterization of relatively compact sets in lp (, b)[END_REF]. We stress that we do not restrict our study to the homogenization method of the bidomain model with nonlinear ionic function of FitzHugh-Nagumo type but also with physiological ionic function of Luo-Rudy type. Moreover, the approach presented herein can be extended to electropermeabilization models. We cite for instance [START_REF] Habib Ammari | Towards monitoring critical microscopic parameters for electropermeabilization[END_REF] where a dynamical homogenization scheme is obtained from a physiological cell model and [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF] where a conductivity dependent macroscopic tissue model is for the first time derived from first principles. Note that thanks to homogenization, the resulting macroscopic bidomain model describes averaged intra and extracellular potential by a nonlinear anisotropic reaction-diffusion system. The cardiac tissue is then considered (at the macroscopic level) as the superposition of two anisotropic continuous media: the intra-and extracellular spaces, coexisting along with the cell membrane, at each point of the tissue. The most substantial mathematical description of the bidomain model is found in the review paper by Henriquez [START_REF] Craig | Simulating the electrical behavior of cardiac tissue using the bidomain model[END_REF], which presents a formal definition of the model from its origins in the core conductor model, and outlines many of the approximations that can be made under certain assumptions. The plan of this paper is outlined as follows. The microscopic problem and the main assumptions used for homogenization are presented in Section 2 and the main result is stated. In Section 3, existence of weak solutions to the microscopic problem is proved based on a Faedo-Galerkin approach, a priori estimates and a compactness argument. In Section 4, some estimates on the solutions of the microscopic problem are obtained and the microscopic problem is formulated using the unfolding operator. The passage to the limit using compactness and the unfolding method are established in Section 5. Then in Section 6, the macroscopic bidomain equations are recuperated from the limit equations obtained in Section 5 and the cell problem is decoupled. Finally, in Section 7, a microscopic bidomain model with physiological ionic model is homogenized to obtain the corresponding macroscopic model.

The microscopic bidomain model

We first list in the following paragraphs the assumptions used in sections 3, 4, 5 and 6.

Assumptions on the domain. For our model we assume that Ω (the cardiac tissue) is a bounded open subset of R 3 with smooth boundary ∂Ω. The cardiac tissue is composed of two connected regions, the intracellular Ω i,ε and the extracellular Ω e,ε . These two regions are separated by an active membrane surface Γ ε = ∂Ω i,ε ∩ ∂Ω e,ε . Here ε > 0 is the small dimensionless parameter which is proportional to the ratio between the micro scale of the length of the cells and the macro scale of the length of the cardiac fibers. Following the standard ap- Therefore, the physical region Ω occupied by the heart is decomposed into the intra-and extracellular domains Ω j,ε for j = i, e that can be simply obtained by intersecting Ω with Y j,k,ε for j = i, e, i.e.:

Ω j,ε = Ω ∩ k∈Z 3 Y j,k,ε . Similarly, Γ ε = Ω ∩ k∈Z 3 Γ k,ε .
One can observe that the domain Ω i,ε may be considered as a perforated domain obtained from Ω by removing the perforations which correspond to the extracellular domain Ω e,ε . The same observation holds for the extracellular domain. The boundary Γ is a smooth manifold such that Γ ε is smooth and connected. Furthermore, Ω j,ε are both assumed to be connected bounded domains in R 3 so that a Poincaré-Wirtinger inequality is satisfied in both domains. (We refer the reader to the geometrical hypothesis H p in [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF] for such domains.)

Assumptions on the data. The electric properties of the tissue are described by the intracellular u i,ε and extracellular u e,ε electric potentials. Herein, u j,ε : Ω j,ε → R for j = i, e, and v ε := (u i,ε -u e,ε ) | Γε : Γ ε → R is known as the transmembrane potential and satisfies a dynamic condition on Γ ε involving the auxiliary function w ε : Γ ε → R (the so called gating variable). The following coupled reaction-diffusion system forms the microscopic bidomain model: for j = i, e (see e.g. [START_REF] Sundnes | Computing the electrical activity in the heart[END_REF][START_REF] Ying | A multilevel adaptive approach for computational cardiology[END_REF]):

-div (M j,ε ∇u j,ε ) = 0 in Ω j,ε,T := (0, T ) × Ω j,ε , (1a)
ε(∂ t v ε + I ion (v ε , w ε ) -I app,ε ) = I m on Γ ε,T := (0, T ) × Γ ε , (1b)
I m = -M i,ε ∇u i,ε • µ i = M e,ε ∇u e,ε • µ e on Γ ε,T , (1c)
∂ t w ε -H(v ε , w ε ) = 0 on Γ ε,T . (1d)
We augment [START_REF] Allaire | Homogenization and two-scale convergence[END_REF] with no-flux boundary conditions

(M j (x)∇u j,ε ) • µ j = 0 on (0, T ) × (∂Ω j,ε \ Γ ε ), j ∈ {e, i}, (2)
and appropriate initial conditions for the transmembrane potential and gating variable

v ε (0, •) = v 0,ε (•), w ε (0, •) = w 0,ε (•) on Γ ε . (3)
The conductivity tensors, the ionic functions, the source term and the initial data satisfy the following assumptions:

(E.1) The conductivity of the tissue is represented by scaled symmetric Lipschitz continuous tensors M i,ε (x) = M i (x, x/ε) and M e,ε (x) = M e (x, x/ε) satisfying (the ellipticity and periodicity conditions): there exists constants m 1 , m 2 > 0 such that for j = i, e

m 1 |ζ| 2 ≤ M j (x, ξ)ζ • ζ ≤ m 2 |ζ| 2 , (4a)
M j (x, ξ + e k ) = M j (x, ξ), (4b)
for all (x, ξ) ∈ Ω × Y j and for all ζ ∈ R 3 . Furthermore, note that µ j are the exterior unit normals to the boundaries of Ω j,ε , for j = i, e respectively, and µ i = -µ e on Γ ε .

(E.

2) The ionic current I ion (u, w) is assumed to be decomposed into I 1,ion (u) and I 2,ion (w), where I ion (u, w) = I 1,ion (u) + I 2,ion (w). Furthermore, the function I 1,ion : R → R is considered as a C 1 function, and the functions I 2,ion : R → R and H : R 2 → R are considered as linear functions. Also, we assume that there exists r ∈ (2, +∞) and constants

α 1 , α 2 , α 3 , L > 0, l ≥ 0 such that 1 α 1 |v| r-1 ≤ |I 1,ion (v)| ≤ α 1 |v| r-1 + 1 ,
and

I 2,ion (w)v -α 2 H(v, w)w ≥ α 3 |w| 2 , (5)
Ĩ1,ion : z → I 1,ion (z) + Lz + l is strictly increasing on R with lim z→0 Ĩ1,ion (z)/z = 0 (6a)

and ∀ z, s ∈ R ( Ĩ1,ion (z) -Ĩ1,ion (s))(z -s) ≥ 1 C (1 + |z| + |s|) r-2 |z -s| 2 . ( 6b
)
Remark 2.1. One can easily show that: I 1,ion (0) = -l, I 1,ion (0) = -L and I 1,ion (z) ≥ -L for all z ∈ R.

Remark 2.2. The function H in the ODE of (1)-( 2) and the function I ion , may correspond to one of the simplified models for the membrane and ionic currents. We mention, for instance, the Mitchell-Schaeffer membrane model [START_REF] Colleen | A two-current model for the dynamics of cardiac membrane[END_REF]

H(v, w) = w ∞ (v/v p ) -w R m c m η ∞ (v/v p ) , (7a)
I ion (v, w) = v p R m v v p η 2 - v 2 (1 -v/v p )w v 2 p η 1 , (7b)
where the dimensionless time constant and state variable constant are respectively given by

η ∞ (s) = η 3 for s < η 5 , η 4 otherwise, w ∞ (s) = 1 for s < η 5 , 0 otherwise.
The quantity R m is the surface resistivity of the membrane, and v p , η 1 , η 2 , η 3 , η 4 , η 5 are given parameters. A simpler choice for the membrane kinetics is given by the widely known FitzHugh-Nagumo model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF], often used by researchers to avoid computational difficulties. In this case,

H(v, w) = av -bw, (8a)
I ion (v, w) = λv(1 -v)(v -θ) + (-λw) := I 1,ion (v) + I 2,ion (w), (8b)
where a, b, λ, θ are given parameters with a, b ≥ 0, λ < 0 and 0 < θ < 1. According to the Mitchell-Shaeffer and FitzHugh-Nagumo models, the most appropriate value is r = 4, which means that the non-linearity I ion is of cubic growth at infinity (recall that in the Mitchell-Shaeffer membrane model, the gating variable w is bounded in L ∞ ). Assumptions (5), ( 6) are automatically satisfied by any cubic polynomial I ion with positive leading coefficient. This is indeed the case for the FitzHugh-Nagumo model but not for the Mitchell-Shaeffer model.

(E.

3) There exists a constant C independent of ε such that the source term I app,ε satisfies the following bound:

ε 1/2 I app,ε L 2 (Γ ε,T ) ≤ C. (9)
Furthermore, I app is the weak limit of the corresponding unfolding sequence.

(E.4) The initial data v 0,ε and w 0,ε satisfy

ε 1/r v 0,ε L r (Γε) + ε 1/2 v 0,ε L 2 (Γε) + ε 1/2 w 0,ε L 2 (Γε) ≤ C, (10)
for some constant C independent of ε. Moreover, v 0,ε and w 0,ε are assumed to be traces of uniformly bounded sequences in C 1 ( Ω).

Observe that the equations in (1) are invariant under the change of u i,ε and u e,ε into u i,ε +k; u e,ε + k, for any k ∈ R. Hence, we may impose the following normalization condition:

Ωe,ε u e,ε (t, x) dx = 0 for a.e. t ∈ (0, T ). [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] Finally, We end this section by stating the main results of the paper as given in the following theorems.

Theorem 2.1 (Microscopic Bidomain Model). Assume conditions (E.1), ..., (E.4) hold.

Then the microscopic bidomain problem (1), ( 2), (3) possesses a unique weak solution in the sense of Definition 3.1.

Theorem 2.2 (Macroscopic Bidomain Model). A sequence of solutions (u i,ε , u e,ε , w ε ) ε of the microscopic system (1)-(3)(obtained in Theorem 2.1) converges to a weak solution (u i , u e , w) 1) (Ω T ) and w ∈ C(0, T ; L 2 (Ω)), of the macroscopic problem

with v = u i -u e , u i , u e ∈ L 2 (0, T ; H 1 (Ω)), v ∈ L 2 (0, T ; H 1 (Ω))∩L r (Ω T ), ∂ t v ∈ L 2 (0, T ; (H 1 (Ω)) )+ L r/(r-
|Γ|∂ t v -div (M i (x)∇u i ) + |Γ|I ion (v, w) = |Γ|I app in Ω T , (12a)
|Γ|∂ t v + div (M e (x)∇u e ) + |Γ|I ion (v, w) = |Γ|I app in Ω T , (12b)
∂ t w -H(v, w) = 0 in Ω T . ( 12c
)
supplemented with no-flux boundary conditions, representing an insulated cardiac tissue

(M j (x)∇u j ) • n = 0 on Σ T := ∂Ω × (0, T ), j ∈ {e, i}, (13)
and appropriate initial conditions in Ω, namely v 0 and w 0 ∈ L 2 (Ω), for the transmembrane potential and gating variable

v(0, x) = v 0 (x), w(0, x) = w 0 (x). ( 14
)
Herein, n is the outward unit normal to the boundary of Ω and the tensors M i and M e are defined by

M j := Y j M j + M j ∇ y f j for j = i, e
, where the components f k,j of f j (k = 1, 2, 3) are the corrector functions, solutions of the cell problems

       -∇ y • (M j ∇ y f k,j ) = -∇ y • (M j e k ) in Y j , M j ∇ y f k,j • µ j = M j e k • µ j on Γ, Y j f k,j = 0, f k,j Y -periodic.

Existence of solutions to the microscopic model

This section is devoted to proving existence of solutions to the microscopic bidomain model for fixed ε > 0. The existence proof is based on the Faedo-Galerkin method, a priori estimates, and the compactness method. We start with a weak formulation of the microscopic model.

Definition 3.1 (Weak Formulation). A solution of problem (1), (2), (3) is a four tuple (u i,ε , u e,ε , v ε , w ε ) such that u i,ε ∈ L 2 (0, T ; H 1 (Ω i,ε )), u e,ε ∈ L 2 (0, T ; H 1 (Ω e,ε )), v ε = (u i,ε - u e,ε ) | Γε ∈ L 2 (0, T ; H 1/2 (Γ ε )) ∩ L r (Γ ε,T ), w ε ∈ L 2 (Γ ε,T )), ∂ t v ε , ∂ t w ε ∈ L 2 (Γ ε,T
), and satisfying the following weak formulation for a.e. t ∈ (0, T )

Γε ε∂ t v ε ϕ ds(x) + j=i,e Ω j,ε M j,ε (x)∇u j,ε • ∇ϕ j dx + Γε εI ion (v ε , w ε )ϕ ds(x) = Γε εI app,ε ϕ ds(x), (15)
Γε ∂ t w ε ζ ds(x) - Γε H(v ε , w ε )ζ ds(x) = 0, (16)
for all ϕ j ∈ H 1 (Ω j,ε ) with

ϕ := (ϕ i -ϕ e ) | Γε ∈ H 1/2 (Γ ε ) ∩ L r (Γ ε ) for j = i, e and ζ ∈ L 2 (Γ ε ).
We prove now Theorem 2.1.

Proof. In this proof, we will remove the ε-dependence in the solution (v ε , u i,ε , u e,ε , w ε ) for simplification of notation. To prove existence of weak solutions, we use a Faedo-Galerkin approach and a priori estimates. For this sake, we first carefully construct an appropriate basis for our systems.

Step 1: Construction of the basis We first consider functions φ ∈ C 0 ( Ωj,ε ) and we define the inner product denoted . Then we let V 0,j denote the completion of C 0 ( Ωj,ε )

under the norm induced by the inner product •, • V 0,j . Similarly, for functions φ, φ ∈ C 1 ( Ωj,ε ), we define the inner product denoted

•, • V 1,j by: Θ, Θ V 1,j := Ω j,ε M j,ε ∇φ • ∇ φdx + Γε φ| Γε φ| Γε ds + Γε ∇ Γε φ • ∇ Γε φds,
where ∇ Γε denotes the tangential gradient operator on Γ ε and we let V 1,j denote the completion of C 1 ( Ωj,ε ) under the norm induced by the inner product •, • V 1,j . We note that the following injections hold:

V 0,j ⊂ L 2 (Ω j,ε ), and V 1,j ⊂ H 1 (Ω j,ε ).
Moreover, the injection from V 1,j into V 0,j is continuous and compact. We refer the reader to [START_REF] Gal | Well-posedness and long time behavior of the non-isothermal viscous cahn-hilliard equation with dynamic boundary conditions[END_REF][START_REF] Racke | The cahn-hilliard equation with dynamic boundary conditions[END_REF] for similar approaches.

It follows from a well-known result (see e.g. [START_REF] Temam | Infinite-dimensional dynamical systems in mechanics and physics[END_REF] p. 54) that the closed bilinear form a(Θ, Θ) := Θ, Θ V 1,j defines a strictly positive self-adjoint unbounded operator

B j : D(B j ) = {Θ ∈ V 1,j : B j Θ ∈ V 0,j } → V 0,j such that, for any Θ ∈ V 1,j , we have B j Θ, Θ V 0,j = a(Θ, Θ). Thus, for k ∈ N, we take a complete system of eigenfunctions Θ k,j = φ k,j ψ k,j k of the problem B j Θ k,j = λ k Θ k,j in V 0,j with Θ k,j ∈ D(B j
), and ψ k,j = φ k,j | Γε where φ k,j and ψ k,j are regular enough. Moreover, the eigenvectors {Θ k,j } k , form an orthogonal basis in V 1,j and V 0,j , and they may be assumed to be normalized in the norm of V 0,j . Since

C 1 ( Ωj,ε ) ⊂ V 1,j ⊂ H 1 (Ω j,ε ), and C 1 ( Ωj,ε ) is dense in H 1 (Ω j,ε ), then V 1,j is dense in H 1 (Ω j,ε ) for the H 1 norm. Therefore, {Θ k,j } k is a basis in H 1 (Ω j,ε ) for the H 1 norm.
On the other hand, we consider a basis {ζ k } k , k ∈ N that is orthonormal in L 2 (Γ ε ) and orthogonal in H 1 (Γ ε ) and we set the spaces

T j,n = span{Θ 1,j , • • • , Θ n,j }, T j,∞ = ∞ n=1 T j,n , K n = span{ζ 1 , • • • , ζ n }, K ∞ = ∞ n=1 K n ,
where T ∞ and K ∞ are dense subspaces of V 1,j and H 1 (Γ ε ) respectively.

Step 2: Construction of approximate solutions For any n ∈ N, we are looking for functions of the form

u j,n ūj,n = n k=1 d j,k (t) φ j,k ψ j,k , j = i, e, with φ j,k | Γε = ψ j,k and w n = n k=1 c k (t)ζ k (x), (17)
solving the approximate regularized problem:

(ε + δ n ) Γε ∂ t ūi,n ψ i ds(x) -ε Γε ∂ t ūe,n ψ i ds(x) + δ n Ω i,ε ∂ t u i,n φ i dx = Γε (-I ion (v n , w n ) + I app,ε )ψ i ds(x) - Ω i,ε M i,ε (x)∇u i,n • ∇φ i dx (18) -ε Γε ∂ t ūi,n ψ e ds(x) + (ε + δ n ) Γε ∂ t ūe,n ψ e ds(x) + δ n Ωe,ε ∂ t u e,n φ e dx = Γε (I ion (v n , w n ) -I app,ε )ψ e | Γε ds(x) - Ωe,ε M e,ε (x)∇u e,n • ∇φ e dx (19) Γε ∂ t w n ζ ds(x) = Γε H(v n , w n )ζ ds(x), (20)
where

δ n = 1 n , Θ j = φ j ψ j ∈ T j,n , for j = i, e and ζ ∈ K n . The terms δ n Γε ∂ t ūj,n ψ j ds(x)
and δ n Ω j,ε ∂ t u j,n φ j dx, j = i, e were added to overcome the degeneracy in [START_REF] Cioranescu | The periodic unfolding method in homogenization[END_REF].

We aim to apply the standard existence theorems for ODEs. For this purpose, if n fixed, we choose Θ i = Θ k,i , Θ e = Θ k,e and ζ = ζ k , 1 ≤ k ≤ n and we substitute the expressions [START_REF] Cioranescu | The periodic unfolding method in perforated domains[END_REF] to the unknowns u i,n , ūi,n , u e,n , ūe,n , and w n . The ODE system, that we obtain, has as unknowns the column vectors

d i = {d i,k } n k=1 , d e = {d e,k } n k=1 and c = {c k } n k=1 .
It can be written as follows:

           (δ n + ε) Āii d i -ε Āie d e + δ n A ii d i = F i (t, d i , d e , c) -ε Āie d i + (δ n + ε) Āee d e + δ n A ee d e = F e (t, d i , d e , c) G c (t) = H(t, d i , d e , c), (21)
where the (k, l) entry of the matrix Āmj , m, j

= i, e is ψ m,k , ψ j,l L 2 (Γε) , for 1 ≤ k, l ≤ n, the (k, l) entry of the matrix Ājj , j = i, e is φ j,k , φ j,l L 2 (Ω j,ε ) , the (k, l) entry of the matrix G is ζ k , ζ l L 2 (Γε)
and where the right hand side vectors F i , F e and H assemble the right hand sides of the equations given in ( 18)- [START_REF] Dobberschütz | Homogenization of a diffusion-reaction system with surface exchange and evolving hypersurface[END_REF]. Note that by the orthonormality of the basis, the matrix

G = ζ k , ζ l L 2 (Γε) 1≤k,l≤n = I n×n ,
is the identity matrix. Furthermore, the first two systems of equations in system ( 21) can be written in the following form:

δ n Āii + A ii 0 0 Āee + A ee + ε Āii -Āie -ĀT ie Āee d i d e = F i F e . ( 22
)
Now making use of the orthonormality of the bases in the spaces V 0,j , the matrices Ājj + A jj , for j = i, e, are equal to the identity n × n matrix I n×n . So system ( 22) may be written as

M d i d e = F i F e , where M = δ n I n×n 0 0 I n×n + ε Āii -Āie -ĀT ie Āee . (23)
In order to write

d i d e = M -1 F i F e
, one needs to prove that the matrix M is invertible.

For this sake, it is enough to prove that the matrix

N := Āii -Āie -ĀT ie Āee is positive semi- definite. Let d = d i d e
, where

d i = (d i,1 , • • • , d i,n ) T ∈ R n and d e = (d e,1 , • • • , d e,n ) T ∈ R n . Then d T Nd = d T i Āii d i -2d T i Āi,e d e + d T e
Āee d e So we have

d T Nd = Γε k,l [d i,k d i,l ψ ik ψ il -2d i,k d e,l ψ ik ψ el + d e,k d e,l ψ ek ψ el ] ds(x) = Γε l d i,l ψ il - l d e,l ψ el 2 ds(x) ≥ 0.
Thus the matrix M is symmetric positive definite, hence invertible. Consequently, the whole system (21) can be written as a system of ordinary differential equations in the form y (t) = f (t, y(t)).

Moreover, the problem that we obtained is supplemented with initial conditions

u i,n (0, x) = u 0,i,n (x) := n l=1 d i,l (0)φ i,l (x), ūi,n (0, x) = ū0,i,n (x) := n l=1 d i,l (0)ψ i,l (x), d i,l (0) := u i,0 ūi,0 , Θ i,l V i,0 , u e,n (0, x) = u 0,e,n (x) := n l=1 d e,l (0)φ e,l (x), ūe,n (0, x) = ū0,e,n (x) := n l=1 d e,l (0)ψ e,l (x), d e,l (0) := u e,0 ūe,0 , Θ e,l V e,0 , w n (0, x) = w 0,n (x) := n l=1 c n,l (0)ζ l (x), c n,l (0) := (w 0 , ζ l ) L 2 (Γε) . (24)
Proceeding exactly as in Ref. [START_REF] Bendahmane | Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue[END_REF], we prove that the entries of F i , F e and H are Carathéodory functions bounded by L 1 functions and we obtain the local existence on the interval [0, t ) of the Faedo-Galerkin solutions u i,n , u e,n , v n and w n .

The global existence of the Faedo-Galerkin solutions is a consequence of the n-independent estimates that are derived in the next section. For more details, consult Ref. [START_REF] Bendahmane | Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue[END_REF].

Step 3: Energy estimates Note that the Galerkin solutions satisfy the following weak formulations:

Γε ε∂ t v n ϕ n ds(x) + i,e Γε δ n ε∂ t ūj,n φj,n ds(x) + i,e Ω j,ε δ n ∂ t u j,n ϕ j,n dx + i,e Ω j,ε M j,ε (x)∇u j,n • ∇ϕ j,n dx + Γε εI ion (v n , w n )ϕ n ds(x) = Γε εI app,ε ϕ n ds(x), (25)
Γε ∂ t w n e n ds(x) - Γε H(v n , w n )e n ds(x) = 0, ( 26
)
where the functions ϕ j,n (t, x) := n l=1 b j,n,l (t)φ j,l (x), e n (t, x) := n l=1 z n,l (t)ζ l (x) and ϕ n := φi,n -φe,n for some given absolutely continuous coefficients b j,n,l (t), z n,l (t) for j = i, e. Herein, φj,n is the trace of ϕ j,n on Γ ε for j = i, e. Now, substituting ϕ j,n = u j,n and e n = εα 2 w n in ( 25) and ( 26), respectively, integrating over (0, s) for s ∈ (0, T ] and summing the resulting equations, one obtains upon using ( 5) and ( 6), Young's inequality, the uniform ellipticity of M j,ε and the L 2 bound on I app,ε :

1 2 ε 1/2 v n (s) 2 L 2 (Γε) + α 2 ε 1/2 w n (s) 2 L 2 (Γε) + i,e ε 1/2 δ 1/2 n ūj,n (s) 2 L 2 (Γε) + i,e δ 1/2 n u j,n (s) 2 L 2 (Ω j,ε ) +m 1 i,e ∇u j,n 2
L 2 (Ω j,ε,s ) + ε Ĩ1,ion (v n )v n L 1 (Γε,s) ≤ 1 2 ε 1/2 v 0,n 2
L 2 (Γε) + w 0,n 2
L 2 (Γε) + i,e ε 1/2 δ 1/2 n ū0,j,n 2 L 2 (Γε) + i,e δ 1/2 n u 0,j,n 2 L 2 (Ω j,ε ) + s 0 Γε εI app,ε v n ds(x) dt - s 0 Γε εI 2,ion (w n ) v n ds(x) dt +α 2 ε s 0 Γε H(v n , w n ) w n ds(x) dt + s 0 Γε ε(L v n + l) v n ds(x) dt ≤ C s 0 ( ε 1/2 v n 2 L 2 (Γε) + ε 1/2 w n 2 L 2 (Γε) ) dt + 1 , (27)
for some constant C independent of n and ε. Note that in the sequel C is a generic constant whose value can change from one line to another. One obtains from [START_REF] Grandelius | The cardiac bidomain model and homogenization[END_REF], the following inequality

ε 1/2 v n (s) 2 L 2 (Γε) + ε 1/2 w n (s) 2 L 2 (Γε) ≤ C s 0 ( ε 1/2 v n 2 L 2 (Γε) + ε 1/2 w n 2 L 2 (Γε) ) dt + 1 .
Hence, by an application of Gronwall's inequality, one gets for a.e. t ∈ (0, T ),

ε 1/2 v n (t) 2 L 2 (Γε) + ε 1/2 w n (t) 2 L 2 (Γε) ≤ C.
Therefore,

ε 1/2 v n L ∞ (0,T ;L 2 (Γε)) + ε 1/2 w n L ∞ (0,T ;L 2 (Γε)) ≤ C.
Exploiting this last inequality along with [START_REF] Grandelius | The cardiac bidomain model and homogenization[END_REF], one obtains

√ εv n L ∞ (0,T ;L 2 (Γε)) + j=i,e δ n √ εū j,n L ∞ (0,T ;L 2 (Γε)) (28)
+

j=i,e

δ n u j,n L ∞ (0,T ;L 2 (Ω j,ε )) + √ εw n L ∞ (0,T ;L 2 (Γε)) ≤ C, ε Ĩ1,ion (v n )v n L 1 (Γ ε,T ) + j=i,e ∇u j,n L 2 (Ω j,ε,T ) ≤ C, (29)
ε 1/r v n L r (Γ ε,T ) ≤ C, (30)
√ εv n L 2 (Γ ε,T ) + √ εw n L 2 (Γ ε,T ) ≤ C, (31)
for some constant C > 0 not depending on n and ε. Moreover, one can obtain some uniform estimates on the time derivatives as follows. Substitute ϕ i,n = ∂ t u i,n and ϕ e,n = ∂ t u e,n in [START_REF] Gal | Well-posedness and long time behavior of the non-isothermal viscous cahn-hilliard equation with dynamic boundary conditions[END_REF], and integrate in time to deduce

ε Γ ε,T |∂ t v n | 2 ds(x) dt + j=i,e Γ ε,T δ n ε(∂ t ūj,n ) 2 ds(x) dt + j=i,e Ω j,ε,T δ n (∂ t u j,n ) 2 dx dt + j=i,e Ω j,ε,T M j,ε (x)∇u j,n • ∇(∂ t u j,n ) dx dt + ε Γ ε,T I 1,ion (v n )∂ t v n ds(x) dt + ε Γ ε,T I 2,ion (w n )∂ t v n ds(x) dt = ε Γ ε,T I app,ε ∂ t v n ds(x) dt. (32)
Now, set

P M j,ε (s) = 1 2 Ω j,ε M j,ε ∇u j,n • ∇u j,n dx and I 1 (s) = s 0 I 1,ion (v)dv. Observe that Ω j,ε,T M j ∇u j,n • ∇(∂ t u j,n ) dx dt = T 0 ∂ t P M j,ε dt = P M j,ε (T ) -P M j,ε (0),
and

Γ ε,T I 1,ion (v n )∂ t v n ds(x) dt = T 0 ∂ t Γε I 1 (v n )ds(x) dt = Γε I 1 (v n (T, y)) ds(x) - Γε I 1 (v n (0, y)) ds(x).
Using this and Young's inequality, one gets from (32)

ε Γ ε,T |∂ t v n | 2 ds(x) dt + j=i,e Γ ε,T δ n ε(∂ t ūj,n ) 2 ds(x) dt + j=i,e Ω j,ε,T δ n (∂ t u j,n ) 2 dx dt + j=i,e P M j,ε (T ) + ε Γε I 1 (v n (T, y)) ds(x) ≤ -ε Γ ε,T I 2,ion (w n )∂ t v n ds(x) dt + P M j,ε (0) + Γε I 1 (v n (0, y)) ds(x) + ε Γ ε,T I app,ε ∂ t v n ds(x) dt ≤ ε 2 Γ ε,T |∂ t v n | 2 ds(x) dt + Cε Γ ε,T |w n | 2 ds(x) dt + P M j,ε (0) + ε Γε I 1 (v n (0, y)) ds(x) + C Γ ε,T |I app,ε | 2 ds(x) dt, (33)
for some constant C > 0 not depending on ε and δ n (recall that 0 < ε ≤ 1).Note that by the monotonicity of Ĩ1,ion (see ( 6)), one can obtain

ε Γε I 1 (v n (T, y)) + L 2 |v n (T, y)| 2 + lv n (T, y) ds(x) ≥ 0. Finally, use j=i,e |P M j,ε (0)| + ε Γε I 1 (v n (0, y)) ds(x) ≤ C j=i,e Ω j,ε |∇u j,n (0, x)| 2 dx + ε Γε |v n (0, y)| r ds(x),
(for some constant C > 0) and estimates ( 28) and ( 29) to get from (33)

ε Γ ε,T |∂ t v n | 2 ds(x) dt + j=i,e Γ ε,T δ n ε(∂ t ūj,n ) 2 ds(x) dt + j=i,e Ω j,ε,T δ n (∂ t u j,n ) 2 dx dt ≤ C, (34)
for some constant C > 0. Hence, one has the estimate

√ ε ∂ t v n L 2 (0,T ;L 2 (Γε)) + i,e δ n √ ε ∂ t ūj,n L 2 (0,T ;L 2 (Γε)) ≤ C, (35)
for some constant C > 0 not depending on n. Also, exploiting the structure of (26) along with estimate (31), one obtains

√ ε∂ t w n L 2 (Γ ε,T ) ≤ C, (36)
for some constant C > 0 independent of n. The above estimates are not sufficient since estimates on the L 2 norms of the intracellular and extracellular potentials are needed in Ω i,ε and Ω e,ε respectively. Due to the compatibility condition [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF], an application of Poincaré's inequality (see for instance [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF]) implies that

u e,n L 2 (0,T ;H 1 (Ωe,ε)) ≤ C. (37)
Furthermore, making use of the trace inequality as stated in [START_REF] Allaire | Two-scale convergence on periodic surfaces and applications[END_REF], one has

ε ūe,n 2 L 2 (Γ ε,T ) ≤ C( u e,n 2
L 2 (Ω e,ε,T ) + ε ∇u e,n 2
L 2 (Ω e,ε,T ) ),
and consequently

ε ūe,n 2 L 2 (Γ ε,T ) ≤ C. ( 38
) Moreover, having ε v n 2 L 2 (Γ ε,T ) ≤ C, there holds ε ūi,n 2 L 2 (Γ ε,T ) ≤ C. (39)
Finally, making use of this last inequality, of ( 29) and of Lemma C.2 in [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF], one gets

u i,n 2
L 2 (Ω i,ε ) ≤ Cε ūi,n 2 L 2 (Γε) + Cε 2 ∇u i,n 2
L 2 (Ω i,ε ) ≤ C. Therefore, the following estimate holds u i,n L 2 (0,T ;H 1 (Ω i,ε )) ≤ C. (40)
The next step is to show that the local solution constructed above can be extended to the whole time interval [0, T ) (independent of n) but this can be done using the above estimates as in Ref. [START_REF] Bendahmane | Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue[END_REF], so we omit the details.

Step 4: Passage to the limit and existence of solutions From ( 40) and [START_REF] Racke | The cahn-hilliard equation with dynamic boundary conditions[END_REF], it is easy to see that v n , ūj,n are bounded in L 2 (0, T ; H 1/2 (Γ ε )) for j = i, e. This is a consequence of the fact that the trace of a function in H 1 is a function in H 1/2 and of the continuity of the trace map. Moreover, we deduce from [START_REF] Luo | A dynamic model of the cardiac ventricular action potential. i. simulations of ionic currents and concentration changes[END_REF], [START_REF] Sundnes | Computing the electrical activity in the heart[END_REF] and [START_REF] Tartar | Quelques remarques sur l'homogénéisation[END_REF] the uniform bound on v n + (-1) j √ δ n ūj,n in L 2 (Γ ε,T ) for j = i, e. Recall that by the Aubin-Lions compactness criterion, the injection

W = {u ∈ L 2 (0, T ; H 1/2 (Γ ε )) and ∂ t u ∈ L 2 (0, T ; H -1/2 (Γ ε ))} ⊂ L 2 (Γ ε,T ))
is compact. Therefore, we can assume there exist limit functions u i,ε , u e,ε , v ε , w ε with v ε = (u i,ε -u e,ε ) | Γε := ūi,ε -ūe,ε on Γ ε,T such that as n → ∞ (for fixed ε and up to an unlabeled subsequence)

                                 v n + (-1) j √ δ n ūj,n → v ε a.e. in Γ ε,T , strongly in L 2 (Γ ε,T ), and weakly in L 2 (0, T ; H 1/2 (Γ ε )) for j = i, e, u j,n u j,ε weakly in L 2 (0, T ; H 1 (Ω j,ε )) for j = i, e, v n → v ε a.e. in Γ ε,T , strongly in L 2 (Γ ε,T ), w n → w ε a.e. in Γ ε,T , strongly in L 2 (Γ ε,T ), I 1,ion (v n ) → I 1,ion (v ε ) a.e. in Γ ε,T and weakly in L r/(r-1) (Γ ε,T ), ∂ t v n ∂ t v ε weakly in L 2 (Γ ε,T ) and δ n ∂ t ūj,n 0 in D (0, T ; L 2 (Γ ε )) for j = i, e, ∂ t w n ∂ t w ε weakly in L 2 (Γ ε,T ), δ n ∂ t u j,n 0 in D (0, T ; L 2 (Ω j,ε )) for j = i, e. (41)
Keeping in mind ( 41), ( 28) and [START_REF] Arsen | Mathematical problems in elasticity and homogenization[END_REF] we infer, by letting n → ∞ in ( 18), ( 19) and [START_REF] Dobberschütz | Homogenization of a diffusion-reaction system with surface exchange and evolving hypersurface[END_REF],

ε Γε ∂ t v ε ϕ ds(x) + i,e Ω j,ε M ε j (x)∇u j,ε • ∇ϕ j dx + ε Γε I ion (v ε , w ε )ϕ ds(x) = ε Γε I app,ε ϕ ds(x), (42)
Γε ∂ t w ε φ ds(x) - Γε H(v ε , w ε )φ ds(x) = 0, (43)
for all

ϕ j ∈ H 1 (Ω j,ε ), j = i, e, with ϕ := (ϕ i -ϕ e ) | Γε ∈ H 1/2 (Γ ε ) ∩ L r (Γ ε ) and φ ∈ L 2 (Γ ε ).
Step 5: Uniqueness. Let (u i,ε,1 , u e,ε,1 , w ε,1 ) and (u i,ε,2 , u e,ε,2 , w ε,2 ) be two weak solutions satisfying ( 42)-( 43), with v ε,k = (u i,ε,k -u e,ε,k ) | Γε for k = 1, 2 and with "data" v ε,0 = v ε,0,1 , w ε,0 = w ε,0,1 and v ε,0 = v ε,0,2 , w ε,0 = w ε,0,2 respectively. Note that the following equations hold for all test functions

ϕ j ∈ L 2 (0, T ; H 1 (Ω j,ε )), j = i, e, with ϕ := (ϕ i -ϕ e ) | Γε ∈ L 2 (0, T ; H 1/2 (Γ ε )) ∩ L r (Γ ε,T ) and φ ∈ L 2 (Γ ε,T ): Γε,t ε∂ t (v ε,1 -v ε,2 )ϕ ds(x) ds + i,e Ω(j,ε,t) M j,ε (x)∇(u j,ε,1 -u j,ε,2 ) • ∇ϕ j dx ds + Γε,t ε(I ion (v ε,1 , w ε,1 ) -I ion (v ε,2 , w ε,2 ))ϕ ds(x) ds = 0, Γε,t ∂ t (w ε,1 -w ε,2 )φ ds(x) ds - Γε,t (H(v ε,1 , w ε,1 ) -H(v ε,2 , w ε,2
))φ ds(x) ds = 0, for 0 < t ≤ T . Substituting ϕ j = (u j,ε,1 -u j,ε,2 ) and φ = w ε,1 -w ε,2 in the two equations above, then adding the resulting ones, we arrive at

1 2 Γε ε |(v ε,1 -v ε,2 )(t)| 2 + |(w ε,1 -w ε,2 )(t)| 2 ds(x) - 1 2 Γε ε|v ε,1,0 -v ε,2,0 | 2 + |w ε,1,0 -w ε,2,0 | 2 ds(x) + j=i,e t 0 Ω j,ε M j,ε ∇(u j,ε,1 -u j,ε,2 ) • ∇(u j,ε,1 -u j,ε,2 ) dx ds + t 0 Γε ε(I ion (v ε,1 , w ε,1 ) -I ion (v ε,2 , w ε,2 ))(v ε,1 -v ε,2 ) ds(x) ds = t 0 Γε (H(v ε,1 , w ε,1 ) -H(v ε,2 , w ε,2 ))(w ε,1 -w ε,2 ) ds(x) dt.
Now using (4), one has for j = i, e

t 0 Ω j,ε M j,ε ∇(u j,ε,1 -u j,ε,2 ) • ∇(u j,ε,1 -u j,ε,2 ) dx ds ≥ 0.
Also, by (6) there holds

t 0 Γε ε(I 1,ion (v ε,1 ) -I 1,ion (v ε,2 ))(v ε,1 -v ε,2 ) ds(x) dt ≥ -Lε t 0 Γε (v ε,1 -v ε,2 ) 2 ds(x) dt.
Moreover, exploiting the linearity of H(v, w) and I 2,ion (w), and using Young's inequality one can deduce

1 2 Γε ε |(v ε,1 -v ε,2 )(t)| 2 + |(w ε,1 -w ε,2 )(t)| 2 ds(x) ≤ C t 0 Γε ε |v ε,1 -v ε,2 | 2 + |w ε,1 -w ε,2 | 2 ds(x) ds + 1 2 Γε ε|v ε,1,0 -v ε,2,0 | 2 + |w ε,1,0 -w ε,2,0 | 2 ds(x),
for some constant C > 0. Finally, an application of Gronwall's inequality yields

Γε ε |(v ε,1 -v ε,2 )(t)| 2 +|(w ε,1 -w ε,2 )(t)| 2 ds(x) ≤ C Γε ε|v ε,1,0 -v ε,2,0 | 2 +|w ε,1,0 -w ε,2,0 | 2 ds(x).
for some constant C > 0. This completes the uniqueness proof.

Convergence of solutions to the macroscopic problem

This section consists in preparing the ground for the passage to the limit as ε → 0. First, some a priori estimates are obtained on the solutions of the microscopic problem. Then, the unfolding operator for perforated domains and the boundary unfolding operator are introduced and some of their properties are recalled. Finally, the microscopic problem is written in an equivalent formulation, the so called "unfolded" formulation, making use of the unfolding operators.

4.1. Energy estimates for the microscopic solutions. The following estimates follow from the estimates on the Faedo-Galerkin solutions obtained in the previous section. Lemma 4.1. Assume that conditions (E.1),...,(E.1) and (1) and (2) hold. Then there exist constants c 1 , c 2 , c 3 , c 4 > 0, not depending on ε such that

√ εv ε L ∞ (0,T ;L 2 (Γε)) + √ εw ε L ∞ (0,T ;L 2 (Γε)) ≤ c 1 , (44)
j=i,e

u j,ε L 2 (0,T ;H 1 (Ω j,ε )) ≤ c 2 (45) ε 1/r v ε L r (Γ ε,T ) ≤ c 3 and ε (r-1)/r I 1,ion (v ε ) L r/(r-1) (Γ ε,T ) ≤ c 4 . ( 46
) If v ε,0 ∈ H 1/2 (Γ ε ) ∩ L r (Γ ε )
, then there exists a constant c 5 > 0 not depending on ε such that In this subsection, we view the domains Ω j,ε , j = i, e as perforated domains and we define the unfolding operator T j ε , j = i, e following the same notation as in [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF]. First, we define the following sets in R 3 (see Figure 2):

√ ε∂ t v ε L 2 (Γ ε,T ) ≤ c 5 . (47)
Ξ ε = { ∈ Z n , ε( + Y ) ⊂ Ω}, Ωε = interior ∈Ξε ε( + Ȳ ) , Ωj ε = interior ∈Ξε ε( + Ȳj ) , j = i, e Ωε,T = (0, T ) × Ωε , Ωj ε,T = (0, T ) × Ωj ε , Γε := {y ∈ Γ ε : y ∈ Ωε }, Λ ε = Ω \ Ωε , Λ ε,T = (0, T ) × Ω \ Ωε .
Secondly, we recall the definition of the time dependent unfolding operator in perforated domains.

Definition 4.1. For any function φ Lebesgue-measurable on (0, T ) × Ω j,ε , the unfolding operator is defined by

T j ε (φ)(t, x, y) = φ t, ε x ε Y + εy a.e. for (t, x, y) ∈ Ωε,T × Y j 0 a.e. for (t, x, y) ∈ (0, T ) × Ω \ Ωε × Y j , (48)
where [•] denotes the Gauβ-bracket. Observe that the function ε[ x ε ] represents the lattice translation point of the ε-cellular medium containing x. For the sake of completeness, we recall some properties of the aforementioned operator and we refer the reader to [START_REF] Cioranescu | The periodic unfolding method in perforated domains[END_REF][START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF] for details.

Proposition 4.2. For p ∈ [1, ∞), the operator T j ε is linear and continuous from L p ((0, T ) × Ω j,ε ) to L p (Ω T × Y j ). For every φ ∈ L 1 ((0, T ) × Ω j,ε ) and v, w ∈ L p ((0, T ) × Ω j,ε ), there holds

(1) T j ε (vw) = T j ε (v)T j ε (w), (2)
Ω×Y j T j ε (φ)(t, x, y) dx dy = Ωj,ε φ(t, x) dx, (3)
T j ε (w) L p (Ω×Y j ) = w1 Ωj ε L p (Ω j,ε ) ≤ w L p (Ω j,ε
) . Furthermore, since the dynamic equations are defined on the surface Γ ε , we resort to the use of the boundary unfolding operator, developed in [START_REF] Cioranescu | The periodic unfolding method in perforated domains[END_REF][START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF] and defined as follows (recall that Γ ε,T = (0, T ) × Γ ε ):

T b ε : L 2 (Γ ε,T ) → L 2 (Ω T × Γ) such that T b ε u(t, x, y) = u(t, ε([ x ε ] + y)), a.e. for (t, x, y) ∈ Ωε,T × Γ, 0,
a.e. for (t, x, y)

∈ (0, T ) × Ω \ Ωε × Γ. ( 49
)
We also list herein some properties of the boundary unfolding operator as given in [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF].

Proposition 4.3. The boundary unfolding operator has the following properties:

(1)

T b ε is a linear operator. (2) T b ε (φψ) = T b ε (φ)T b ε (ψ), ∀φ, ψ ∈ L p (Γ ε,T ), p ∈ (1, +∞). (3) For every φ ∈ L 1 (Γ ε,T ), we have the following integration formula Γε φ(t, x) ds(x) = 1 ε Ω×Γ T b ε (φ)(t,
x, y) dx ds(y).

(4) For every φ ∈ L p ( Γε,T ) with p ∈ (1, +∞), one has

T b ε (φ) L p (Ω T ×Γ) = ε 1/p φ L p ((0,T )× Γε) .
(5) For every ϕ ∈ D(Ω T × Γ) and w ∈ W 1,1 (0, T ; L 1 (Γ ε )), the following integration by parts formula holds

T 0 Ω×Γ T b ε (∂ t w)T b ε (ϕ) ds(y) dx dt = - T 0 Ω×Γ T b ε (w)T b ε (∂ t ϕ) ds(y) dx dt.
Remark 4.1. Note that the last property (which is not listed in [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF]) is a direct consequence of the integration by parts formula:

T 0 Γε ∂ t wϕ ds(x) dt = - T 0 Γε w∂ t ϕ ds(x) dt,
and the integration formula in property (3) of Proposition 4.3. Now, in order to make use of the unfolding method in the homogenization of the microscopic problem, we rewrite the corresponding equations ( 42) and [START_REF] Vogt | A homogenization theorem leading to a Volterra-integrodifferential equation for permeation chromatography[END_REF] in the "unfolded" form. We have the following identities:

Ω j,ε,T M j,ε ∇u j,ε • ∇ϕ j dx dt = Ω T ×Y j T j ε (M j,ε )T j ε (∇u j,ε )T j ε (∇ϕ j ) dx dy dt + (0,T )×Λ j ε M j,ε ∇u j,ε • ∇ϕ j dx dy dt := Ω T ×Y j T j ε (M j,ε )T j ε (∇u j,ε )T j ε (∇ϕ j ) dx dy dt + r 1 , Γ ε,T εI app,ε ϕ ds(x) dt = Γε,T εI app,ε ϕ ds(x) dt + Γ ε,T ∩Λ ε,T εI app,ε ϕ ds(x) dt = Ω T ×Γ T b ε (I app,ε )T b ε (ϕ) ds(y) dx dt + Γ ε,T ∩Λ ε,T εI app,ε ϕ ds(x) dt := Ω T ×Γ T b ε (I app,ε )T b ε (ϕ) ds(y) dx dt + r 2 , Γ ε,T ε∂ t v ε ϕ ds(x) dt = Γε,T ε∂ t v ε ϕ ds(x) dt + Γ ε,T ∩Λ ε,T ε∂ t v ε ϕ ds(x) dt := Ω T ×Γ T b ε (∂ t v ε )T b ε (ϕ) ds(y) dx dt + r 3 , Γ ε,T εI 1,ion (v ε )ϕ ds(x) dt = Γε,T εI 1,ion (v ε )ϕ ds(x) dt + Γ ε,T ∩Λ ε,T εI 1,ion (v ε )ϕ ds(x) dt := Ω T ×Γ T b ε I 1,ion (v ε ) T b ε (ϕ) ds(y) dx dt + r 4 , Γ ε,T εI 2,ion (w ε )ϕ ds(x) dt = Γε,T εI 2,ion (w ε )ϕ ds(x) dt + Γ ε,T ∩Λ ε,T εI 2,ion (w ε )ϕ ds(x) dt := Ω T ×Γ T b ε (I 2,ion (w ε ))T b ε (ϕ) ds(y) dx dt + r 5 = Ω T ×Γ I 2,ion (T b ε (w ε ))T b ε (ϕ) ds(y) dx dt + r 5
Due to the above equalities, one obtains the following equivalent "unfolded" formulation of (42):

Ω T ×Γ T b ε (∂ t v ε )T b ε (ϕ) ds(y) dx dt + i,e Ω T ×Y j T j ε (M j,ε )T j ε (∇u j,ε )T j ε (∇ϕ j ) dx dy dt -λ Ω T ×Γ T b ε (w ε )T b ε (ϕ) ds(y) dx dt + Ω T ×Γ T b ε I 1,ion (v ε ) T b ε (ϕ) ds(y) dx dt = Ω T ×Γ T b ε (I app,ε )T b ε (ϕ) ds(y) dx dt + r 2 -r 5 -r 4 -r 3 -r 1 .
(50) Similarly, the "unfolded" formulation of ( 43) is given by:

Ω T ×Γ T b ε (∂ t w ε )T b ε (φ) ds(y) dx dt - Ω T ×Γ H(T b ε (v ε ), T b ε (w ε ))T b ε (φ) ds(y) dx dt = -ε Γ ε,T ∩Λ ε,T ∂ t w ε φ ds(x) dt + ε Γ ε,T ∩Λ ε,T H(v ε , w ε )φ ds(x) dt := r 6 + r 7 .
(51)

"Unfolding" compactness

In this section, we establish the passage to the limit in (50) and (51). First, note that by estimates ( 44)-(47) obtained above one has

r 1 , • • • , r 7 → 0 as ε → 0.
Also, by regularity of the test functions ϕ and φ, there holds

T b ε ϕ → ϕ and T b ε φ → φ strongly in L 2 (Ω T × Γ),
and T j ε ϕ j → ϕ j strongly in L 2 (Ω T × Y j ). Now, consider Ψ j and θ j in D(Ω T ) and ψ j = ψ j (ξ) in H 1 per (Y j ) such that Y j ψ j = 0 and test equation (50) with functions ϕ ε j = Ψ j + εθ j ψ j,ε where ψ j,ε (x) = ψ( x ε ) (see for e.g. [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF]). Since ∇ϕ ε j = ∇ x Ψ j + εψ j,ε ∇ x θ j + θ j (∇ ξ ψ j,ε ), and thanks to Proposition 2.8 in [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF] (see also [START_REF] Graf | A convergence result for the periodic unfolding method related to fast diffusion on manifolds[END_REF]), there holds

T j ε (ϕ ε j ) → Ψ j strongly in L 2 (Ω T × Y j ), T j ε (θ j ψ j,ε ) → θ j (x)ψ j (ξ) strongly in L 2 (Ω T × Y j ), T j ε (∇ϕ ε j ) → ∇Ψ j + θ j ∇ ξ ψ j strongly in L 2 (Ω T × Y j ), and T b ε (ϕ ε ) → Ψ strongly in L 2 (Ω T × Γ), (52)
where

ϕ ε = (ϕ ε i -ϕ ε e )| Γ ε,T and Ψ = (Ψ i -Ψ e )| Ω T ×Γ .
Hence, to establish the passage to the limit in (50) and (51), we need to verify that the remaining terms of the equations are weakly convergent. Now, making use of estimate (29), there exist u j ∈ L 2 (0, T ; H 1 (Ω)) and ûj ∈ L 2 (0, T ; L 2 (Ω, H 1 per (Y j ))) such that, up to a subsequence (see for instance theorem 3.12 in [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF]), the following hold

T j ε (u ε j ) u j weakly in L 2 (0, T ; L 2 (Ω, H 1 (Y j ))), T j ε (∇u ε j ) ∇u j + ∇ ξ ûj weakly in L 2 (Ω T × Y j ).
Thus, since T j ε (M j,ε ) → M j a.e. in Ω × Y j , one obtains (recall the strong convergence (52))

i,e Ω T ×Y j T j ε (M j,ε )T j ε (∇u ε j )T j ε (∇ϕ j ) dy dx dt → i,e Ω T ×Y j M j (∇u j + ∇ ξ ûj )(∇Ψ j + θ j ∇ ξ ψ j (ξ)) dy dx dt as ε → 0.

Furthermore, since

T b ε (w ε ) L 2 (Ω T ×Γ) ≤ |Y | 1/2 ε 1/2 w ε L 2 (Γ ε,T ) ≤ C, then up to a subsequence T b ε w ε w in L 2 (Ω T × Γ).
Consequently, by linearity of I 2,ion ,

Ω T ×Γ I 2,ion (T b ε (w ε ))T b ε (ϕ) ds(y) dx dt → Ω T ×Γ I 2,ion (w)Ψ ds(y) dx dt.
Similarly, exploiting assumption (9), one obtains

Ω T ×Γ T b ε (I app,ε )T b ε (ϕ) ds(y) dx dt → Ω T ×Γ I app Ψ ds(y) dx dt.
In order to establish the convergence of T b ε (∂ t v ε ), first note that

T b ε (∂ t v ε ) L 2 (Ω T ×Γ) ≤ |Y | 1/2 ε 1/2 ∂ t v ε L 2 Γ ε,T ≤ C. So there exists h ∈ L 2 (Ω T × Γ) such that T b ε (∂ t v ε ) h weakly in L 2 (Ω T × Γ)
. By a classical argument, one can identify h to ∂ t v Therefore,

Ω T ×Γ T b ε (∂ t v ε )T b ε (ϕ ε ) ds(y) dx dt → Ω T ×Γ ∂ t vΨ ds(y) dx dt.
It remains to obtain the passage to the limit in the term containing the ionic function I 1,ion . Indeed due to the nonlinearity, it is difficult to pass to the limit in I 1,ion on the microscopic membrane surface and one needs to establish the passage to the limit in:

lim ε→0 T 0 Ω Γ I 1,ion (T b ε v ε )T b ε (ϕ ε )ds(y)dxdt.
By regularity of ϕ ε , we have

T b ε ϕ ε → Ψ strongly in L r ((0, T ) × Ω × Γ). It remains to show the weak convergence of I 1,ion (T b ε v ε ) to I 1,ion (v) in L r/(r-1) (Ω T ×Γ).
Therefore, we show the strong convergence of

T b ε v ε to v in L 2 (Ω T × Γ).
Then, by the properties of I 1,ion we actually obtain the strong convergence of I 1,ion (T b ε v ε ) to I 1,ion (v) in L q (Ω T × Γ) for all q ∈ [1, r/(r -1)). For this sake, we make use of Kolmogorov-Riesz-type compactness criterion for the space L p (Ω, B) that can be found as Corollary 2.5 in [START_REF] Gahn | A characterization of relatively compact sets in lp (, b)[END_REF].

Proposition 5.1. Let Ω ⊂ R n be an open and bounded set. Let p ∈ [1, ∞), B be a Banach space and F ⊂ L p (Ω, B). Then F is relatively compact in L p (Ω, B) iff (i) for every measurable set C ⊂ Ω, the set { C f dx : f ∈ F } is relatively compact in B, (ii) for all δ > 0 and z ∈ R n and z i ≥ 0, i = 1, . . . , n, there holds

sup f ∈F τ zσ f -f L p (Ω z δ ,B) → 0, for z → 0,
where Ω z δ := {x ∈ Ω δ : x + z ∈ Ω δ } and Ω δ := {x ∈ Ω : dist(x, ∂Ω) > δ}, (iii) for δ > 0, there holds sup f ∈F Ω\Ω δ |f (x)| p dx → 0 for δ → 0.

Fisrt, we prove an estimate on the space translates of the transmembrane potential v ε that is needed later to obtain an estimate on the space translate of T b ε (v ε ). Now, we fix open sets K and K such that K ⊂⊂ K ⊂⊂ Ω, and we let z ∈ R with |z| < dist(K , ∂Ω).

We have the following lemma Lemma 5.2. Let l ∈ Z 3 and ε > 0 such that ε|l| ≤ |z|. Then the following estimate holds:

ε v ε (t, x + εl) -v ε (t, x) 2 L 2 (Γ ε,K ) ≤ Cεl, (53)
where Γ ε,K = Γ ε ∩ K and C is a positive constant.

For simplicity of notation, we use

τ εl v(t, x) := v(t, x + εl).
Proof. In this proof, we consider ϕ j ∈ H 1 (Ω j,ε ) with suppϕ j ⊂ K, for j = i, e. We use the translations of ϕ j , j = i, e i.e. ϕ j (x -εl) as test functions in the variational formulation [START_REF] Cioranescu | The periodic unfolding method in homogenization[END_REF]. In the resulting equation, we make the substitution x → x + εl, and we exploit the periodicity of the domain to get

ε Γε∩K ∂ t (τ εl v ε )(ϕ i -ϕ e ) ds(x) + j=i,e Ω j,ε ∩K τ εl M j,ε ∇τ εl u j,ε • ∇ϕ j dx (54) +ε Γε∩K I ion (τ εl v ε , τ εl w ε )(ϕ i -ϕ e ) ds(x) = Γε∩K τ εl I app,ε ϕ ds(x).
Noting that the last equality is valid for test functions with support in K , let η ∈ D(K ) be a cutoff function for K, with 0 ≤ η ≤ 1, η = 1 in K and zero outside K . We test the variational equation for (τ εl u j,ε -u j,ε ), j = i, e with

ϕ j = η 2 (τ εl u j,ε -u j,ε ), j = i, e,
we get

ε 2 d dt Γε∩K η 2 (τ εl v ε -v ε ) 2 ds(x) (55) + j=i,e Ω j,ε ∩K τ εl M j,ε ∇τ εl u j,ε -M j,ε ∇u j,ε • ∇ η 2 (τ εl u j,ε -u j,ε ) dx +ε Γε∩K η 2 I ion (τ εl v ε , τ εl w ε ) -I ion (v ε , w ε ) (τ εl v ε -v ε ) ds(x) = Γε∩K (τ εl I app,ε -I app,ε )η 2 (τ εl v ε -v ε ) ds(x).
First, we break up the second term in (55) as follows:

j=i,e Ω j,ε ∩K

τ εl M j,ε ∇τ εl u j,ε -M j,ε ∇u j,ε • ∇ η 2 (τ εl u j,ε -u j,ε ) dx = j=i,e Ω j,ε ∩K η 2 τ εl M j,ε ∇τ εl u j,ε -M j,ε ∇u j,ε • ∇ τ εl u j,ε -u j,ε dx + j=i,e Ω j,ε ∩K 2η τ εl M j,ε ∇τ εl u j,ε -M j,ε ∇u j,ε • τ εl u j,ε -u j,ε ∇ηdx := T 1 + T 2 , (56)
and we estimate T 1 exploiting the ellipticity of M j,ε given in (4):

T 1 = j=i,e Ω j,ε ∩K η 2 M j,ε ∇τ εl u j,ε -∇u j,ε • ∇ τ εl u j,ε -u j,ε dx + j=i,e Ω j,ε ∩K η 2 τ εl M j,ε -M j,ε ∇(τ εl u j,ε ) • ∇τ εl u j,ε -∇u j,ε dx ≥ 0 -C j=i,e η 2 (τ εl M j,ε -M j,ε ) L ∞ (Ω j,ε ∩K ) ∇u j,ε 2 L 2 (Ω j,ε ∩K ) ≥ -ε|l|C j=i,e ∇u j,ε 2
L 2 (Ω j,ε ∩K ) . (57)
In the last inequality, we used the mean value theorem to obtain:

τ εl M j,ε -M j,ε L ∞ (Ω j,ε ∩K ) ≤ ε|l| j=i,e ∇M ε j L ∞ (Ω j,ε ∩K ) ≤ ε|l|C,
for some constant C > 0. Moreover, by regularity of η, Cauchy-Schwarz and boundedness of M j,ε , we get the following estimate on T 2 :

|T 2 | ≤ Cε|l| j=i,e u j,ε H 1 (Ω j,ε ∩K ) , (58)
for some constant C > 0. On the other hand, the third term of (55) may be divided into two terms by making use of (8b) as follows:

ε Γε∩K η 2 I ion (τ εl v ε , τ εl w ε ) -I ion (v ε , w ε ) (τ εl v ε -v ε ) ds = ε Γε∩K η 2 I 1,ion (τ εl v ε ) -I 1,ion (v ε ) (τ εl v ε -v ε ) ds +ε Γε∩K η 2 I 2,ion (τ εl w ε ) -I 2,ion (w ε ) (τ εl v ε -v ε ) ds := T 3 + T 4 . (59)
By monotonicity ( 6), we estimate T 3 :

T 3 ≥ -εL η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) . (60)
In addition, using the definition of η and the linearity of I 2,ion (5), Cauchy-Schwarz and Young's inequalities, T 4 can be estimated by:

|T 4 | ≤ εC( τ εl w ε -w ε 2 L 2 (Γε∩K ) + η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) ), (61)
for some constant C > 0. Furthermore, the source term in (55) satisfies the following inequality:

Γε∩K (τ εl I app,ε -I app,ε )η 2 (τ εl v ε -v ε ) ds(x) ≤ Cε I app,ε L 2 (Γε) η(τ εl v ε -v ε ) L 2 (Γε∩K ) . ( 62
)
Gathering all these estimates, one obtains

ε d dt η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) ≤ C 1 ε|l| + C 2 ε η(τ εl w ε -w ε ) 2 L 2 (Γε∩K ) + η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) . (63)
By a similar argument, one can also obtain from ( 16),

ε d dt η(τ εl w ε -w ε ) 2 L 2 (Γε∩K ) ≤ C 3 ε η(τ εl w ε -w ε ) 2 L 2 (Γε∩K ) + η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) . ( 64
)
By Grönwall's inequality applied to the sum of ( 63) and (64), we obtain

ε η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) + ε η(τ εl w ε -w ε ) 2 L 2 (Γε∩K ) ≤ e C 4 t C 1 ε|l|t + ε η(τ εl v 0,ε -v 0,ε ) 2 L 2 (Γε∩K ) + ε η(τ εl w 0,ε -w 0,ε ) 2 L 2 (Γε∩K ) ≤ C(T ) ε|l| + ε τ εl v 0,ε -v 0,ε 2
L 2 (Γε∩K ) + ε τ εl w 0,ε -w 0,ε 2 L 2 (Γε∩K )
. Now using the assumption on v 0,ε and w 0,ε , one obtains

ε η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) + ε η(τ εl w ε -w ε ) 2 L 2 (Γε∩K ) ≤ Cε|l| Furthermore, noting that η(τ εl v ε -v ε ) 2 L 2 (Γε∩K ) + η(τ εl w ε -w ε ) 2 L 2 (Γε∩K ) ≥ τ εl v ε -v ε 2 L 2 (Γε∩K) + τ εl w ε -w ε 2 L 2 (Γε∩K)
, one can conclude that (53) holds. Now, we state and prove the strong convergence of T b ε (v ε ) to v. Lemma 5.3. The following convergence holds:

T b ε (v ε ) → v strongly in L 2 (Ω T × Γ), as ε → 0. Moreover, I 1,ion (T b ε v ε ) → I 1,ion (v) strongly in L q (Ω T × Γ), for q ∈ [1, r/(r -1)), as ε → 0
Proof. The proof of the lemma is similar to the proof of Theorem 14 in [START_REF] Gahn | Homogenization of reaction-diffusion processes in a two-component porous medium with nonlinear flux conditions at the interface[END_REF] but herein the domain is not the union of scaled and translated reference cells. The proof is based on Proposition 5.1. Condition (iii) follows from estimate (46), since

Ω\Ω δ |T b ε (v ε )| 2 dx ≤ |Ω \ Ω δ | r-2 r Ω |T b ε (v ε )| r dx 2 r ≤ C|Ω \ Ω δ | r-2 r .
To prove condition (i), consider a measurable set A ⊂ Ω, and define

v ε A (t, y) = A T b ε (v ε )(t, x, y) dx, for a.e. t ∈ (0, T ), y ∈ Γ.
The a priori estimates obtained on v ε imply that v ε A is bounded in L 2 ((0, T ), H 1/2 (Γ)) ∩ H 1 ((0, T ), L 2 (Γ)). Then by Aubin-Lions Lemma, the sequence is relatively compact in L 2 ((0, T ), L 2 (Γ)). Using the properties of the unfolding operator (see Proposition 4.3) and estimate [START_REF] Ying | A multilevel adaptive approach for computational cardiology[END_REF] one can easily find a constant C > 0 such that

T 0 v ε A 2 L 2 (Γ) ≤ C.
By a similar argument and making use of the estimate on ε 1/2 ∂ t v ε , one can also show that

∂ t v ε A L 2 (Γ T ) ≤ C, for some positive constant C.
On the other hand, to obtain a uniform estimate on the L 2 (0, T ;

H 1/2 (Γ)), we first observe that v ε A 2 H 1/2 (Γ) = v ε A 2 L 2 (Γ) + |v ε A | 2 H 1/2 0 (Γ)
.

Based on the previous estimates, we only need to bound the H 1/2 0 seminorm and this is done as follows. First, we have by Cauchy-Schwarz and Fubini

|v ε A | 2 H 1/2 0 (Γ) ≤ |A| Ωε Γ Γ |v ε (t, ε[ x ε ] + εy 1 ) -v ε (t, ε[ x ε ] + εy 2 )| 2 |y 1 -y 2 | 3
ds(y 1 ) ds(y 2 )dx.

We note that this is equivalent to writing:

|v ε A | 2 H 1/2 0 (Γ) ≤ |A| Ωε |v ε (t, ε[ x ε ] + ε•)| 2 H 1/2 0 (Γ)
dx.

Since v ε = (u i,ε -u e,ε )

Γ and using the triangle inequality, we get

|v ε A | 2 H 1/2 0 (Γ) ≤ 2|A| j=i,e Ωε |u j,ε (t, ε[ x ε ] + ε•)| 2 H 1/2 0 (Γ)
dx.

Now, by the trace inequality which can be found in [START_REF] Allaire | Two-scale convergence on periodic surfaces and applications[END_REF], we find a constant C > 0 such that

|v ε A | 2 H 1/2 0 (Γ) ≤ C|A| j=i,e Ωε u j,ε (t, ε[ x ε ] + ε•) 2 L 2 (Y j ) + ∇ y (u j,ε (t, ε[ x ε ] + ε•)) 2 L 2 (Y j ) dx.
By the chain rule, we have

∇ y u j,ε (t, ε[ x ε ] + εy) = ε∇u j,ε t, ε[ x ε ] + εy . So Ωε ∇ y (u j,ε (t, ε[ x ε ] + εy)) 2 L 2 (Y j ) dx = Ωε ε∇u j,ε (t, ε[ x ε ] + εy) 2 L 2 (Y j ) dx, or equivalently Ωε ∇ y (u j,ε (t, ε[ x ε ] + εy)) 2 L 2 (Y j ) dx = ε 2 Ωε Y j ∇u j,ε (t, ε[ x ε ] + εy) 2 dydx. Now, using Proposition 4.2-(2), we get Ωε ∇ y (u j,ε (t, ε[ x ε ] + εy)) 2 L 2 (Y j ) dx = ε 2 |Y | Ωj,ε ∇u j,ε 2 dx = ε 2 |Y | ∇u j,ε 2
L 2 ( Ωj,ε ) .
One more time, we make use of Proposition 4.2-(2) to obtain:

Ωε Y j u j,ε t, ε[ x ε ] + εy dy dx = |Y | Ωj,ε u j,ε (t, x) dx, and
|v ε A | 2 H 1/2 0 (Γ) ≤ C j=i,e u j L 2 ( Ωj,ε ) + ε 2 ∇u j,ε 2
L 2 ( Ωj,ε ) .
Finally, integrating over (0, T ) and using the a priori estimates (29) on u j,ε , we obtain the required result. It remains to prove condition (ii) of Proposition 5.1 as follows.

Fix ε > 0 and let I ⊂ Z 3 , be an index set such that

Ωε = i∈I ε(Y + i).
Obviously, we have

x ∈ ε(Y + i) ⇔ [ x ε ] = i. For every i ∈ I we divide the cell ε(Y + i) into subsets ε(Y + i) k with k ∈ {0, 1} 3 , defined as follows ε(Y + i) k := x ∈ ε(Y + i) : ε x + { ξ ε }ε ε = ε(i + k) ,
for a given ξ ∈ R 3 such that ξ is O(ε). Then we have the following identity:

ε(Y + i) = k∈{0,1} 3 ε(Y + i) k . Now, we compute τ ξ T b ε (v ε ) -T b ε (v ε ) 2 L 2 ((0,T )×Ω ξ δ ×Γ) = τ ξ T b ε (v ε ) -T b ε (v ε ) 2 L 2 ((0,T )×(Ω ξ δ ∩ Ωε)×Γ) + τ ξ T b ε (v ε ) -T b ε (v ε ) 2 L 2 ((0,T )×(Ω ξ δ ∩ Ωc ε )×Γ) ≤ E 1,ξ,ε + E 2,ξ,ε , where E 1,ξ,ε := τ ξ T b ε (v ε ) -T b ε (v ε ) 2 L 2 ((0,T )× Ωε×Γ)
, and

E 2,ξ,ε := τ ξ T b ε (v ε ) -T b ε (v ε ) 2 L 2 ((0,T )×(Ω ξ δ \ Ωε)×Γ) = τ ξ T b ε (v ε ) 2 L 2 ((0,T )×(Ω ξ δ \ Ωε)×Γ)
.

We first estimate E 1,ξ,ε , making use of the fact that Ωε = i∈I ε(Y + i), and proceeding in a similar way to [START_REF] Dobberschütz | Homogenization of a diffusion-reaction system with surface exchange and evolving hypersurface[END_REF][START_REF] Neuss | Effective transmission conditions for reaction-diffusion processes in domains separated by an interface[END_REF] as follows:

E 1,ξ,ε = i∈I T 0 ε(Y +i) Γ v ε t, ε x + ξ ε + εy -v ε t, ε x ε + εy 2 ds(y) dx dt = i∈I k∈{0,1} 3 T 0 ε(Y +i) k Γ v ε t, ε i + k + ξ ε + εy -v ε (t, εi + εy) 2 ds(y) dx dt ≤ i∈I k∈{0,1} 3 T 0 ε(Y +i) Γ v ε t, ε i + k + ξ ε + εy -v ε (t, εi + εy) 2 ds(y) dx dt ≤ k∈{0,1} 3 T 0 Ωε Γ T b ε v ε t, x + ε k + ξ ε , y -T b ε v ε (t, x , y) 2 ds(y) dx dt
≤ k∈{0,1} 3 ε T 0 Γε v ε t, x + ε ξ ε + k -v ε (t, x) 2 ds(x) dt,
where in the last inequality we used the identity i = x ε and the integration formula of Proposition 4.3-(3). Moreover, using estimate (53), we obtain

E 1,ξ,ε ≤ C(|ξ| + ε).
Therefore, one can conclude that E 1,ξ,ε → 0 as ξ → 0 uniformly in ε, as in [START_REF] Gahn | Homogenization of reaction-diffusion processes in a two-component porous medium with nonlinear flux conditions at the interface[END_REF]. Indeed, to prove that ∀ρ > 0, ∃µ > 0 such that ∀ε > 0, ∀ξ, |ξ| ≤ µ ⇒ E 1 < ρ, one identifies two cases:

(

1) ε < ρ 2C : take µ = ρ 2C , then for ξ < µ, E 1 < ρ. (2) ρ 2C < ε: since ε -1 ∈ N and 1 ε < 2C ρ ≤ [ 2C ρ ] + 1
, there are finitely many values ε such that ε > ρ 2C , and for each such

ε i i = 1, • • • , m, ∃µ i such that ∀ξ, |ξ| < µ i ⇒ E 1 < ρ by continuity of translation in L 2 . Take µ 0 = min{µ, µ 1 , • • • , µ m }.
Then the estimate follows.

Consider now E 2,ξ,ε , and note that

E 2,ξ,ε ≤ τ ξ T b ε (v ε ) 2 L 2 ((0,T )×(Ω δ \ Ωε)×Γ)
.

Observe that, for ε small enough, say ε < ε 0 , Ω δ ⊂ Ωε , so E 2,ξ,ε = 0. On the other hand, for ε 0 < ε < 1, since ε -1 ∈ N, there exist finitely many ε ∈ (ε 0 , 1), say

{ε j } m j=1 , m ∈ N, m < ∞.
Moreover, by continuity of the translation of L 2 functions, for each ρ > 0 there exists for every j, a β(ε j ) such that

E 2,ξ,ε < ρ, ∀|ξ| < β(ε j ). Let β = min{β(ε 1 ), • • • , β(ε m )}, then for all ρ > 0, |ξ| < β ⇒ E 2,ξ,ε < ρ.
Hence, E 2,ξ,ε → 0 as ξ → 0, uniformly in ε. This ends the proof of (ii) in Proposition 5.1.

The following result is therefore obtained:

T b ε (v ε ) → v strongly in L 2 (Ω T × Γ),
as ε → 0. Finally, to prove the convergence of the nonlinear term in the ionic function, first note that from the structure of I 1,ion and using the properties of the boundary unfolding operator, there holds

T b ε (I 1,ion (v ε )) = I 1,ion (T b ε (v ε )), then using the estimate ε (r-1)/r I 1,ion (v ε ) L r/(r-1) (Γ ε,T ) ≤ C, one obtains T b ε (I 1,ion (v ε )) L r/(r-1) (Ω T ×Γ) ≤ |Y | (r-1)/r ε (r-1)/r I 1,ion (v ε ) L r/(r-1)(Γ ε,T ) ≤ C.
Hence, since up to a subsequence

T b ε (v ε ) → v a.e. in Ω T × Γ,
one gets, using the continuity of I 1,ion and a classical result (see Lemma 1.3 in [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]),

I 1,ion (T b ε (v ε )) I 1,ion (v) weakly in L r/(r-1) (Ω T × Γ).
Moreover, using Vitali's theorem, one has the strong convergence of

I 1,ion (T b ε (v ε )) to I 1,ion (v) in L q (Ω T × Γ) for q ∈ [1, r/(r -1)).
Collecting all the convergence results stated above, one obtains the following limiting problem:

|Γ| Ω T ∂ t vΨ dx dt + i,e Ω T ×Y j M j [∇u j + ∇ ξ ûj ][∇Ψ j + θ j ∇ ξ ψ j ] +|Γ| Ω T I 2,ion (w)Ψ dx dt + |Γ| Ω T I 1,ion (v)Ψ dx dt = |Γ| Ω T I app Ψ dx dt. (65)
Similarly, one can easily show that the limit of (51) as ε tends to 0, is given by

|Γ| Ω T ∂ t wφ dx dt -|Γ| Ω T H(v, w)φ dx dt = 0. (66)
6. Macroscopic bidomain model (proof of Theorem 2.2)

The next step is to obtain the weak formulation of the bidomain equations and the cell problem. So one needs to formulate the limit problem in terms of u i and u e alone and hence find an expression of ûi and ûe in terms of u i , u e respectively. First, to determine the cell problem, set in (65) Ψ i , Ψ e and θ e to 0, to get

Ω T ×Y i M i [∇u i + ∇ y ûi ][θ i ∇ y ψ i ] dy dx dt = 0,
which corresponds to the classical cell problem obtained in section 2 and it can be shown that the function ûi can be written in terms of u i as follows (û i is defined up to an additive function in x, see for instance [START_REF] Cioranescu | The periodic unfolding method in domains with holes[END_REF]):

ûi (t, x, y) = f i (t, x, y) • ∇ x u i + f 0,i (t, x) = 3 k=1 ∂u i ∂x k f k,i (t, x, y) + f 0,i (t, x), (67)
where the corrector functions (i.e. the components of the function

f i ) f k,i ∈ L ∞ (Ω T ; H 1 per (Y i )), k = 1, 2, 3, are for a.e. (t, x) ∈ Ω T the solutions of the cell problems        -∇ y • (M i ∇ y f k,i ) = -∇ y • (M i e k ) in Y i , M i ∇ y f k,i • µ i = M i e k • µ i on Γ, Y i f k,i = 0, f k,i Y -periodic. (68)
The existence and uniqueness of the correctors follow by classical arguments from Lax-Milgram theorem (see for instance the remark on p. 13-14 of [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF] or [START_REF] Arsen | Mathematical problems in elasticity and homogenization[END_REF]). Finally, inserting formula (67) into (65) and setting θ i ,Ψ e and θ e to 0, one obtains the weak formulation of the macroscopic bidomain model

|Γ| Ω T ∂ t vΨ i dx dt + Ω T M i ∇u i • ∇Ψ i + |Γ| Ω T I 2,ion (w)Ψ i dx dt +|Γ| Ω T I 1,ion (v)Ψ i dx dt = |Γ| Ω T I app Ψ i dx dt, (69)
where M i is elliptic and defined by

M i := Y i M i + M i ∇ y f i .
Similarly, one can decouple the cell problem in the extracellular domain and define the homogenized conductivity matrix M e .

Remark 6.1. Since the convergence obtained herein is shown up to a subsequence, it is required to prove uniqueness of the macroscopic problem to guarantee the convergence of the whole sequence. Indeed, uniqueness of the macroscopic bidomain model has been obtained for several ionic models, we refer for instance to [START_REF] Bendahmane | Analysis of a class of degenerate reaction-diffusion systems and the bidomain model of cardiac tissue[END_REF][START_REF] Bourgault | Existence and uniqueness of the solution for the bidomain model used in cardiac electrophysiology[END_REF] for the case of phenomenological models of FitzHugh-Nagumo type and to [START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF] for physiological models of Luo-Rudy type.

Unfolding homogenization to physiological models

In this section, we extend the homogenization results obtained in the previous sections to physiological ionic models. So the ordinary differential equation (1d) is replaced by a system of ODEs for the gating variables w l , l = 1, • • • , k and the concentration variable z. The kinetics of a general physiological ionic model may be represented by the functions R, G and I ion that satisfy assumptions (A.1)-(A.3), stated below. It can be verified that those assumptions are satisfied by several gating and ionic concentration variables in Beeler-Reuter or Luo-Rudy ionic models [START_REF] Veneroni | Reaction-diffusion systems for the macroscopic bidomain model of the cardiac electric field[END_REF][START_REF] Go | Reconstruction of the action potential of ventricular myocardial fibres[END_REF][START_REF] Luo | A dynamic model of the cardiac ventricular action potential. i. simulations of ionic currents and concentration changes[END_REF].

R l (v, w) = α l (v)(1 -w l ) -β l (v)w l ( 70
)
where α l and β l , l = 1, • • • , k are positive rational functions of exponentials in v such that:

0 < α l (v), β l (v) ≤ C α,β (1 + |v|). ( 71
) (A.
2) The function I ion : R × R k × (0, +∞) → R has the general form:

I ion (v, w, z) = k l=1 I l ion (v, w l ) + I z ion (v, w, z, ln z) (72)
where

I l ion ∈ C 0 (R × R k ) ∩ Lip(R × [0, 1] k
) and satisfies the condition:

|I l ion (v, w l )| ≤ C 1,I (1 + |w l | + |v|), (73)
and I z ion is such that:

I z ion ∈ C 1 (R × R k × R + × R) ∩ Lip(R × [0, 1] k × R + × R), I z ion (v, w, z, ln z) ≤ C 2,I (1 + |v| + |w| + |z| + ln z), (74)
I z ion (v, w, z, ln z) ≥ C 3,I k l=1 (|v| + w l + w l ln z), (75)
0 < Θ(w) ≤ ∂ ∂ζ I z ion (v, w, z, ζ) ≤ Θ(w), (76)
∂ ∂v I z ion (v, w, z, ζ) ≤ L(w), (77)
∂ ∂w l I z ion ≤ C 4,I (1 + |v| + | ln z|), ∀l = 1, • • • , k, (78) 0
≤ ∂ ∂z I z ion ≤ C 5,I , (79)
where Θ, Θ, L belong to C 0 (R, R + ) and C 1,I , . . . , C 5,I are positive constants.

(A.3) The function G ∈ Lip(R × [0, 1] k × R +
) is given by:

G(v, w, z) = a 1 (a 2 -z) -a 3 I z ion (v, w, z, ln z), (80)
where a 1 , a 2 , a 3 are positive physiological constants that vary from one ion to another. In our case, we only consider z to correspond to the intracellular calcium concentration. Under those assumptions, the microscopic system that we consider is given by:

-div (M j,ε ∇u j,ε ) = 0 in Ω j,ε,T := (0, T ) × Ω j,ε , (81a)

ε(∂ t v ε + I ion (v ε , w ε , z ε ) -I app,ε ) = I m on Γ ε,T := (0, T ) × Γ ε , (81b)
I m = -M i,ε ∇u i,ε • µ i = M e,ε ∇u e,ε • µ e on Γ ε,T , (81c)
∂ t w ε -R(v ε , w ε ) = 0 on Γ ε,T , (81d)
∂ t z ε -G(v ε , w ε , z ε ) = 0 on Γ ε,T . (81e)
We augment (81) with no-flux boundary conditions

(M j,ε (x)∇u j,ε ) • µ j = 0 on (0, T ) × (∂Ω j,ε \ Γ ε ), j ∈ {e, i}, (82)
and appropriate initial conditions for the transmembrane potential, the gating variables and the concentration variable

v ε (0, •) = v 0,ε (•), w ε (0, •) = w 0,ε (•), z ε (0, •) = z 0,ε (•) on Γ ε , (83)
where

v 0,ε ∈ H 1/2 (Γ ε ), z 0,ε ∈ L 2 (Γ ε ) and w 0,ε ∈ L 2 (Γ ε ) k with z 0,ε > c 0 > 0 for some c 0 > 0 and 0 ≤ w l,0,ε ≤ 1 for l = 1, • • • , k.
Analogously to the miscroscopic model with more general FitzHugh-Nagumo dynamics, one has the following existence result.

Theorem 7.1. Suppose that assumptions (A.1)-

(A.3) hold. If v 0,ε ∈ H 1/2 (Γ ε ), z 0,ε ∈ L 2 (Γ ε )
and w 0,ε ∈ L 2 (Γ ε ) k with z 0,ε > c 0 > 0 for some c 0 > 0 and 0 ≤ w l,0,ε ≤ 1 for l = 1, • • • , k, then the microscopic problem (81),(82),(83) possesses a weak solution defined as follows: u i,ε ∈ L 2 (0, T ; H 1 (Ω i,ε )), u e,ε ∈ L 2 (0, T ; H 1 (Ω e,ε )), with Ωe,ε∩Ω u e,ε = 0,

v ε = (u i,ε -u e,ε ) | Γε ∈ L 2 (Γ ε,T ), w ε ∈ (L 2 (Γ ε,T )) k , z ε ∈ L 2 (Γ ε,T ), ∂ t v ε , ∂ t z ε ∈ L 2 (Γ ε,T ), and ∂ t w ε ∈ (L 2 (Γ ε,T ) k such that Γ ε,T ε∂ t v ε ϕ ds(x) dt + j=i,e Ω j,ε,T M j,ε (x)∇u j,ε • ∇ϕ j dx dt + Γ ε,T εI ion (v ε , w ε , z ε )ϕ ds(x) dt = Γ ε,T εI app,ε ϕ ds(x) dt, (84)
Γ ε,T ∂ t w l,ε φ ds(x) dt - Γ ε,T R l (v ε , w ε )φ ds(x) dt = 0, ( 85
)
for l = 1, • • • , k and Γ ε,T ∂ t z ε φ ds(x) dt - Γ ε,T G(v ε , w ε , z ε )φ ds(x) dt = 0, ( 86
)
for all ϕ j ∈ L 2 (0, T ; H 1 (Ω j,ε )) with ϕ := (ϕ i -ϕ e ) | Γε ∈ L 2 (0, T ; H 1/2 (Γ ε )) for j = i, e and φ ∈ L 2 (Γ ε,T ).

The proof of the theorem follows closely the steps done in the case above of more general FitzHugh-Nagumo ionic function type. Using approximation systems and applying a Faedo-Galerkin method in space, one can obtain the existence of a weak solution for the approximation systems (similarly to section 4) then by a passage to the limit, the existence for the microscopic problem is obtained based on some technical results and a series of a priori estimates that are listed in the sequel but their detailed proofs are available in [START_REF] Bendahmane | Mathematical analysis of cardiac electromechanics with physiological ionic model[END_REF]. We also refer to [START_REF] Veneroni | Reaction-diffusion systems for the microscopic cellular model of the cardiac electric field[END_REF] where a fixed point approach was used. First, the recovery variables are shown to satisfy the physiological bounds.

Lemma 7.2. Let w l,ε ∈ C([0, T ], L 2 (Γ ε )) and v ε ∈ H 1 (0, T, L 2 (Γ ε )) such that for all ω ∈ L 2 (Γ ε ): Γε ∂ t w l,ε ω = Γε R l (v ε , w l,ε )ω, ( 87
)
where R(v, w) is defined by (70). Assume that 0 ≤ w l,0,ε ≤ 1 a.e. in Γ ε , then

0 ≤ w l,ε ≤ 1, a.e. in Γ ε,T . (88)
Secondly, one has to make sure that the concentration variable stays positive.

Lemma 7.3. Let z ε ∈ C([0, T ], L 2 (Γ ε )), v ε ∈ H 1 (0, T, L 2 (Γ ε )) and w ε ∈ C([0, T ], L 2 (Γ ε ) k ) such that for all ω ∈ L 2 (Γ ε ): Γε ∂ t z ε ω = Γε G(v ε , w ε , z ε )ω, ( 89
)
where G(v, w, z) satifies assumption (A.6) above. Let z 0 : Ω → (0, +∞) such that:

z 0 ∈ L 2 (Γ ε ), z 0 > 0, a.e. in Γ ε .
Then for a.e.

(t, x) ∈ [0, T ] × Γ ε , z > 0.
Thirdly, the concentration variable and its logarithm ln z ε are proved to be controlled by the norm of v ε in the following sense. Lemma 7.4. Under the same assumptions as Lemma 7.3, the concentration variable z ε satisfies the following estimates for a.e. x ∈ Γ ε , t ∈ (0, T ):

|z ε (t, x)| ≤ C(1 + |z 0,ε (x)| + v ε (x) L 2 (0,t) ), ∀t ∈ [0, T ], (90)
| ln z ε (t, x)| ≤ C(1 + |z 0,ε (x)| + |v ε (t, x)| + v ε (x) L 2 (0,t) ) (91) t 0 |∂ s z ε | 2 ≤ C 1 + |z 0,ε ln z 0,ε | + |z 0,ε | 2 + v ε 2 L 2 (0,t) , (92)
t 0 |ln z ε | 2 ≤ C 1 + |z 0,ε ln z 0,ε | + |z 0,ε | 2 + v ε 2 L 2 (0,t) , (93)
Using the above estimates on z ε and w ε , one can control the L 2 norm of I ion by the L 2 norm of v ε and this result will be later used to reach a uniform in ε estimate on v ε . Lemma 7.5. Under the same conditions of Lemma 7.4, there exists a constant C > 0 (dependent on T ) such that

I ion (v ε , w ε , z ε ) 2 L 2 (Γ ε,T ) ≤ C(1 + v ε 2 L 2 (Γ ε,T ) ). (94)
Based on the previous Lemmata, and proceeding in a similar way as in Section 5, one can easily obtain the following estimates on the solutions to the microscopic problem that are required for the passage to the limit as ε → 0 (the detailed derivation can be found in [START_REF] Bendahmane | Mathematical analysis of cardiac electromechanics with physiological ionic model[END_REF]). Lemma 7.6. There exist constants C 1 , C 2 and C 3 independent of ε such that

max t∈[0,T ] √ εv ε 2 L 2 (Γε) + √ εw ε 2 L 2 (Γε) + √ εz ε 2 L 2 (Γε) ≤ C 1 , (95)
j=i,e

u j,ε L 2 (0,T ;H 1 (Ω j,ε )) ≤ C 2 , (96)
√ ε∂ t v ε L 2 (Γ ε,T ) + √ ε∂ t w ε L 2 (Γ ε,T ) k + √ ε∂ t z ε L 2 (Γ ε,T ) ≤ C 3 , (97)
In order to exploit the unfolding method, the weak formulation is written in its "unfolded" form as in section 4.2 above. Equation (84) becomes:

Ω T ×Γ T b ε (∂ t v ε )T b ε (ϕ) ds(y) dx dt + i,e Ω T ×Y j T j ε (M ε j )T j ε (∇u j,ε )T j ε (∇ϕ j ) dy dx dt + Ω T ×Γ T b ε I ion (v ε , w ε , z ε ) T b ε (ϕ) ds(y) dx dt = Ω T ×Γ T b ε (I app,ε )T b ε (ϕ) ds(y) dx dt + r 10 , (98
) where r 10 is considered as a remainder term which involves integrals over the region Λ ε whose measure tends to zero as ε → 0. Similarly, the "unfolded" formulations of (85) and (86) are given by:

Ω T ×Γ T b ε (∂ t w l,ε )T b ε (φ) ds(y) dx dt - Ω T ×Γ T b ε (R l (v ε , w ε ))T b ε (φ) ds(y) dx dt = r 11 , (99)
for l = 1, • • • , k and Ω T ×Γ T b ε (∂ t z ε )T b ε (φ) ds(y) dx dt - Ω T ×Γ T b ε (G(v ε , w ε , z ε ))T b ε (φ) ds(y) dx dt = r 12 , (100)
where r 11 and r 12 are remainder terms that tend to zero as ε → 0. Now, making use of Lemma 7.6, one can repeat the arguments in section 5 to show that there exist u j ∈ L 2 (0, T ; H 1 (Ω)) and ûj ∈ L 2 (0, T ; L 2 (Ω, H 1 per (Y j ))) such that, up to a subsequence, the following hold

T j ε (u ε j ) u j weakly in L 2 (0, T ; L 2 (Ω, H 1 (Y j ))), T j ε (∇u ε j ) ∇u j + ∇ y ûj weakly in L 2 (Ω T × Y j ).
Thus, one obtains

i,e Ω T ×Y j T j ε (M ε j )T j ε (∇u ε j )T j ε (∇ϕ j ) → i,e Ω T ×Y j M j (∇u j + ∇ y ûj )(∇Ψ j + θ j ∇ y ψ j (y)) as ε → 0.
Furthermore, one can also show that

T b ε (v ε ) L 2 (Ω T ×Γ) + T b ε (w ε ) L 2 (Ω T ×Γ) k + T b ε (z ε ) L 2 (Ω T ×Γ) ≤ C, then up to a subsequence T b ε v ε v in L 2 (Ω T × Γ), T b ε w ε w in L 2 (Ω T × Γ) k , T b ε z ε z in L 2 (Ω T × Γ).
Also, note that due to the a priori estimates on the time derivatives (Lemma 7.6), there exists a constant C > 0 such that

T b ε (∂ t v ε ) L 2 (Ω T ×Γ) + T b ε (∂ t w ε ) L 2 (Ω T ×Γ) k + T b ε (∂ t z ε ) L 2 (Ω T ×Γ) ≤ C, consequently one can show as in section 5, that T b ε (∂ t v ε ) ∂ t v, in L 2 (Ω T × Γ), T b ε (∂ t w ε ) ∂ t w in L 2 (Ω T × Γ) k , T b ε (∂ t z ε ) ∂ t z in L 2 (Ω T × Γ).
Similarly, exploiting assumption (9) on the source term I app,ε , one obtains

Ω T ×Γ T b ε (I app,ε )T b ε (ϕ) ds(y) dx dt → Ω T ×Γ I app Ψ ds(y) dx dt.
It remains to establish the passage to the limit in the nonlinear terms involving the ionic function I ion and the functions R and G appearing in the ODE system. Indeed, making use of assumptions (A.1)-(A.3), of Lemma 7.5 and of Lemma 7.6, there exists a constant C > 0 such that

ε 1/2 I ion (v ε , w ε , z ε ) L 2 (Γ ε,T ) + ε 1/2 R(v ε , w ε ) L 2 (Γ ε,T ) k + ε 1/2 G(v ε , w ε , z ε ) L 2 (Γ ε,T ) ≤ C.
Consequently,

T b ε I ion (v ε , w ε , z ε ) L 2 (Ω T ×Γ) + T b ε R(v ε , w ε ) L 2 (Ω T ×Γ) k + T b ε G(v ε , w ε , z ε ) L 2 (Ω T ×Γ) ≤ C.
Moreover, based on definition 49 of the boundary unfolding operator, one can do the following identifications for a.e. (t, x, y) ∈ Ωε,T × Γ:

T b ε I ion (v ε , w ε , z ε ) = I ion T b ε (v ε ), T b ε (w ε ), T b ε (z ε ) , T b ε R(v ε , w ε ) = R(T b ε (v ε ), T b ε (w ε ) , and T b ε G(v ε , w ε , z ε ) = G T b ε (v ε ), T b ε (w ε ), T b ε (z ε ) .
Hence, there exist functions Ĩion , R and G such that up to a subsequence, the following convergences hold

I ion T b ε (v ε ), T b ε (w ε ), T b ε (z ε ) Ĩion , in L 2 (Ω T × Γ), R T b ε (v ε ), T b ε (w ε ) R, in L 2 (Ω T × Γ) k , and G T b ε (v ε ), T b ε (w ε ), T b ε (z ε ) G, in L 2 (Ω T × Γ).
Therefore, to end the passage to the limit, it remains to relate the functions Ĩion , R, and G to I ion (v, w, z), R(v, w) and G(v, w, z) where v, w and z are the respective limits of T b ε (v ε ), T b ε (w ε ) and T b ε (z ε ). This is done in the following proposition.

Remark 7.1. One possibility is to proceed analogously to Section 6 and prove the strong convergence of T b ε (v ε ), T b ε (w ε ) and T b ε (z ε ). This can be done exactly as in Section 6 for T b ε (v ε ). However, it seems out of reach to prove the strong convergence of T b ε (w ε ) and T b ε (z ε ) by a similar argument. Proposition 7.7. Suppose that assumptions (A.1)-(A.3) are satisfied and let v ε , u j,ε , j = i, e, w ε and z ε be weak solutions of the microscopic system (84)-(86) as given in Theorem 7.1. Then there holds

   Ĩion = I ion (v, w, z), R = R(v, w) G = G(v, w, z),
where v, w, z, Ĩion , R and G are the limits of

T b ε (v ε ), T b ε (w ε ), T b ε (z ε ), T b ε I ion (v ε , w ε , z ε ) , T b ε R(v ε , w ε ) and T b ε G(v ε , w ε , z ε ) respectively.
Proof. Due to assumptions (A.1)-(A.3), in particular the Lipschitz conditions, one can prove that there exists K I > 0 such that

(I ion (v 1 , w 1 , z 1 ) -I ion (v 2 , w 2 , z 2 ))(v 1 -v 2 ) -(R(v 1 , w 1 ) -R(v 2 , w 2 )) • (w 1 -w 2 ) -((G(v 1 , w 1 , z 1 ) -G(v 2 , w 2 , z 2 ))(z 1 -z 2 ) ≥ -K I (|v 1 -v 2 | 2 + |w 1 -w 2 | 2 + |z 1 -z 2 | 2 ).
(101) To obtain the result, we proceed as in [START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Collin | Mathematical analysis and 2-scale convergence of a heterogeneous microscopic bidomain model[END_REF] for 2-scale convergence in nonlinear terms. Using the formulation of the unfolded equations ( 98)-(100) with test functions e -λs u j,ε , e -λs w ε and e -λs z ε respectively, then integrating by parts in time and adding the resulting equations one has

1 2 e -λt T b ε (v ε ) 2 L 2 + 1 2 e -λt k l=1 T b ε (w l,ε ) 2 L 2 + 1 2 e -λt T b ε (z ε ) 2 L 2 - 1 2 e -λt T b ε (v ε,0 ) 2 L 2 - 1 2 e -λt k l=1 T b ε (w l,ε,0 ) 2 L 2 - 1 2 e -λt T b ε (z ε,0 ) 2 L 2 + i,e t 0 e -λs Ω×Y j T j ε (M j,ε )T j ε (∇u j,ε ) • T j ε (∇u j,ε ) dy dx ds + t 0 e -λs [ Ω×Γ I ion (T b ε (v ε ), T b ε (w ε ), T b ε (z ε ))T b ε (v ε ) ds(y) dx + λ 2 T b ε (v ε ) 2 L 2 ]ds + t 0 e -λs [ Ω×Γ -R(T b ε (v ε ), T b ε (w ε )) • T b ε (w ε ) ds(y) dx + λ 2 k l=1 T b ε (w l,ε ) 2 L 2 ] ds + t 0 e -λs [ Ω×Γ -G(T b ε (v ε ), T b ε (w ε ), T b ε (z ε ))T b ε (z ε ) ds(y) dx + λ 2 T b ε (z ε ) 2
A ε -2D ε + E ε + I ε + R ε + G ε ≥ 0. ( 104
)
Now, we set for any positive scalar τ , the following test functions ψ l,ε (t, x) = ψ l,0 (t, x, x ε ) + τ ψ l (t, x, x ε ), θ ε (t, x) = θ 0 (t, x, x ε ) + τ θ(t, x, x ε ), ϕ ε j (t, x) = ϕ 0 j (t, x) + εϕ 1 j (t, x, x ε ) + τ ϕ j (t, x), ϕ ε = (ϕ ε i -ϕ ε e )| Γε . Note that the following convergence results hold strongly in L 2 (Ω T × Γ):

T b ε (ψ l,0 ) → ψ l,0 , T b ε (ψ l ) → ψ l , T b ε (θ 0 ) → θ 0 , T b ε (θ) → θ, T b ε (ϕ ε ) → ϕ 0 + τ ϕ. Moreover T j
ε ϕ ε j → ϕ 0 j + τ ϕ j , T j ε (ϕ 1 j ) → ϕ 1 j (t, x, y), T j ε (∇ϕ ε j ) → ∇(ϕ 0 j + τ ϕ j ) + ∇ y ϕ 1 j , strongly in L 2 (Ω T × Y j ), j = i, e. We pass to the limit in (104), showing the limit of each term separately. Similarly, the limits of R ε and G ε can be obtained to get from inequality (104)

A 0 -2D 0 + E 0 + I 0 + R 0 + G 0 ≥ 0. ( 105
)
This last inequality, being true for any test functions ϕ 0 j , ϕ 1 j , ψ 0 , θ 0 , can be shown to be true by a density argument for u j , ûj , j = i, e, v, w and z. Consequently, one can simplify (105) Finally, repeating the argument of section 6 one can easily decouple the limit equations to get the equations of the macroscopic bidomain model (as (69)):

|Γ|

Ω T ∂ t vΨ i dx dt + Ω T M i ∇u i • ∇Ψ i + |Γ| Ω T I ion (v, w, z)Ψ i dx dt = |Γ| Ω T I app Ψ i dx dt, (109)
where M i is elliptic and defined by

M i := Y i M i + M i ∇ y f i ,
in addition to the corresponding cell problem given by (68). Email address: rtalhouk@ul.edu.lb

Figure 1 .

 1 Figure 1. Left: A 2D section of the simplified periodic network of cells. Right: A 2D section of the reference cell Y .

Figure 2 .

 2 Figure 2. The set Ωe ε in dark red, the set Ωi ε in dark blue and the region Λ ε in dark green and light blue.

4. 2 .

 2 Unfolded formulation of the microscopic problem.

(A. 1 )

 1 Define the function R as R(v, w) := R 1 (v, w 1 ), ..., R k (v, w k ) where R l : R 2 → R are globally Lipschitz continuous functions given by

  (I app,ε )T b ε (v ε ) ds(y) dx ds + t 0 e -λs (r 10 + r 11 + r 12 ) ds.

2 e 2 + 1 2 e 2 L 2 + 2 ke 2 L 2 +- 2 L 2

 2222222222 ), observe that one can take λ large enough so that the following inequality holds1 -λT T b ε (v ε ) -T b ε (ϕ ε ) 2 l,ε ) -T b ε (ψ l,ε ) 2 L -λT T b ε (z ε ) -T b ε (θ ε ) j,ε )T j ε (∇u j,ε -∇ϕ j,ε ) • T j ε (∇u j,ε -∇ϕ j,ε ) dy dx dt + (T b ε (v ε ), T b ε (w ε ), T b ε (z ε )) -I ion (T b ε (ϕ ε ), T b ε (ψ ε ), T b ε (θ ε )) (T b ε (v ε ) -T b ε (ϕ ε )) ds(y) dx + λ 2 T b ε (v ε ) -T b ε (ϕ ε ) 2 b ε (v ε ), T b ε (w ε )) -R(T b ε (ϕ ε ), T b ε (ψ ε ))) • (T b ε (w ε ) -T b ε (ψ ε )) ds(y) dx + λ l=1 T b ε (w l,ε ) -T b ε (ψ l,ε ) 2 L 2 dt + T 0 e -λt Ω×Γ -(G(T b ε (v ε ), T b ε (w ε ), T b ε (z ε )) -G(T b ε (ϕ ε ), T b ε (ψ ε ), T b ε (θ ε ))) (T b ε (z ε ) -T b ε (θ ε )) ds(y) dx + λ 2 T b ε (z ε ) -T b ε (θ ε ) 2 L 2 dt ≥ 0 (103)We want to use (102) to simplify the previous inequality. We introduce the following notation:app,ε )T b ε (v ε ) ds(y) dx dt + -λt (r 10 + r 11 + r 12 )dt, ε )T b ε (ϕ ε ) + T b ε (w ε ) • T b ε (ψ ε ) + T b ε (z ε )T b ε (θ ε ) ds(y) dx + j,ε )T j ε (∇(u j,ε )) • T j ε (∇ϕ j,ε ) dy dx dt, E ε := 1 2 e -λT T b ε (ϕ ε ) 2 L 2 + k l=1 T b ε (ψ l,ε ) 2 L 2 + T b ε (θε) j,ε )T j ε (∇ϕ j,ε ) • T j ε (∇ϕ j,ε ) dy dx dt I ε := I ion (T b ε (v ε ), T b ε (w ε ), T b ε (z ε ))T b ε (ϕ ε ) -I ion (T b ε (ϕ ε ), T b ε (ψ ε ), T b ε (θ ε ))(T b ε (v ε ) -T b ε (ϕ ε )) -λT b ε (v ε )T b ε (ϕ ε ) ds(y) dx + λ 2 T b ε (ϕ ε ) ε ), T b ε (w ε )) • T b ε (ψ ε ) + R(T b ε (ϕ ε ), T b ε (ψ ε )) • (T b ε (w ε ) -T b ε (ψ ε ) -λT b ε (w ε ) • T b ε (ψ ε ) ds(y) dx + λ 2 k l=1 T b ε (ψ l,ε ) 2 L 2 dt, b ε (v ε ), T b ε (w ε ), T b ε (z ε ))T b ε (θ ε ) + G(T b ε (ϕ ε ), T b ε (ψ ε ), T b ε (θ ε )) (T b ε (z ε ) -T b ε (θ ε )) -T b ε (z ε )T b ε (θ ε ) ds(y) dx + λ 2 T b ε (θ ε ) 2 L 2 dt. Substituting (102) into (103), we obtain

A 0 : 2 L 2 (Ω) + z 0 2 L 2 (M 2 L 2 Mϕ 0 + τ ϕ 2 L 2

 022222222 = lim ε→0 A ε 0 + τ ϕ) + w • (ψ 0 + τ ψ) + z(θ 0 + τ θ) dx+ j (∇u j + ∇ y ûj )(∇(ϕ 0 j + τ ∇ϕ) + ∇ y ϕ 1 j ) dy dx, E 0 := lim ε→0 E ε = |Γ| 2 e -λT (ϕ 0 + τ ϕ)(T ) 2 L 2 (Ω) + k l=1 (ψ l,0 + τ ψ l )(T ) 2 L 2 (Ω) + (θ 0 + τ θ)(T ) j (∇(ϕ 0 j + τ ϕ j ) + ∇ y ϕ 1 j ) • (∇(ϕ 0 j + τ ϕ j ) + ∇ y ϕ 1 j ) dy dx, I 0 := lim ε→0 I ε ϕ 0 + τ ϕ) -I ion (ϕ 0 + τ ϕ, ψ 0 + τ ψ, θ 0 + τ θ)(v -ϕ 0 -τ ϕ) ds(y) dx -λ|Γ| Ω v(ϕ 0 + τ ϕ) dx + λ|Γ| 2 (Ω)

I∂

  ion (v + τ ϕ, w + τ ψ, z + τ θ) -Ĩ)ϕ + (R(v + τ ϕ, w + τ ψ) -R) • ψ +(G(v + τ ϕ, w + τ ψ, z + τ θ) -G)θ dxdydt + O(τ 2 ) ≥ 0Dividing by τ and then letting τ tends to 0, we find that for all test functions ϕ, ψ and θ, there holdsT 0 e -λt Ω×Γ I ion (v, w, z) -Ĩ)ϕ + (R(v, w) -R) • ψ +(G(v, w, z) -G)θ dxdydt ≥ 0,Using -ϕ, -ψ and -θ for ϕ, ψ and θ one also getsT 0 e -λt Ω×Γ I ion (v, w, z) -Ĩ)ϕ + (R(v, w) -R) • ψ +(G(v, w, z) -G)θ dxdydt ≤ 0,which gives the result of the proposition.Collecting all the convergence results stated above, one obtains the following limiting problem:|Γ| Ω T ∂ t vΨ dx dt + i,e Ω T ×Y j M j [∇u j + ∇ y ûj ][∇Ψ j + θ j ∇ y ψ j ] +|Γ| Ω T I ion (v,w, z)Ψ dx dt = |Γ| Ω T I app Ψ dx dt, t zφ dx dt -|Γ| Ω T G(v, w, z)φ dx dt = 0. (108)

(

  Mostafa Bendahmane) Institut de mathématiques de Bordeaux (IMB) et l'institut de rythmologie et modélisation cardiaque (Liryc), université de Bordeaux et INRIA-Carmen Bordeaux Sud-Ouest Email address: mostafa.bendahmane@u-bordeaux.fr (Fatima Mroue) Laboratoire de Mathématiques, Ecole Doctorale des Sciences et Technologies, Université Libanaise, Hadat, Département de mathématiques, Centrale Nantes Email address: fatima.mroue@ec-nantes.fr (Mazen Saad) Département de mathématiques, Centrale Nantes and Laboratoire de mathématiques Jean Leray Email address: mazen.saad@ec-nantes.fr (Raafat Talhouk) Laboratoire de Mathématiques, Ecole Doctorale des Sciences et Technologies, Université Libanaise, Hadat

Acknowledgment

We would like to thank the anonymous referees for their constructive and thorough reports. Their suggestions significantly improved our original manuscript.