The cardiac bidomain model and homogenization

We provide a rather simple proof of a homogenization result for the bidomain model of cardiac electrophysiology. Departing from a microscopic cellular model, we apply the theory of two-scale convergence to derive the bidomain model. To allow for some relevant nonlinear membrane models, we make essential use of the boundary unfolding operator. There are several complications preventing the application of standard homogenization results, including the degenerate temporal structure of the bidomain equations and a nonlinear dynamic boundary condition on an oscillating surface.


Introduction
The bidomain model [39,13,38] is widely used as a quantitative description of the electric activity in cardiac tissue. The relevant unknowns are the intracellular (u i ) and extracellular (u e ) potentials, along with the so-called transmembrane potential (v := u i − u e ). In this model, the intra-and extracellular spaces are considered as two separate homogeneous domains superimposed on the cardiac domain. The two domains are separated by the cell membrane creating a discontinuity surface for the cardiac potential. Conduction of electrical signals in cardiac tissue relies on the flow of ions through channels in the cell membrane. In the bidomain model, the celebrated Hodgkin-Huxley [23] framework is used to dynamically couple the intra-and extracellular potentials through voltage gated ionic channels.
The bidomain model can be viewed as a PDE system consisting of two degenerate reaction-diffusion equations involving the unknowns u i , u e , v and two conductivity tensors σ i , σ e . These equations are supplemented by a nonlinear ODE system for the dynamics of the ion channels. The bidomain model is often derived heuristically by interpreting σ i , σ e as some sort of "average" conductivities, applying Ohm's electrical conduction law and the continuity equation (conservation of electrical charge) to the intracellular and extracellular domains [13,38].
Starting from a more accurate microscopic (cell-level) model of cardiac tissue, with the heterogeneity of the underlying cellular geometry represented in great detail, it is possible to heuristically derive the bidomain model (tissue-level) using the multiple scales method of homogenization. This derivation was first carried out in [31]. It should be noted that the microscopic model is in general too complex to allow for full organ simulations, although there have been some very recent efforts in that direction [40]. The complexity of cell-level models, which themselves can be heuristically derived from the Poisson-Nernst-Planck equations [35], motivates the search for simpler homogenized (macroscopic) models. The work [31] assumes, as we do herein, that cardiac tissue can be viewed as a uniformly oriented periodic assembly of cells (see also [14,22]). There have been some attempts to remove this assumption. We refer to [26,25,36] for extensions to somewhat more realistic tissue geometries.
Despite the widespread use of the bidomain model, there are few mathematical rigorous derivations of the model from a microscopic description of cardiac tissue. From a mathematical point of view, rigorous homogenization is often linked to the study of the asymptotic behavior (convergence) of solutions to PDEs with oscillating coefficients. In the literature several approaches have been developed to handle this type of problem, like Tartar's method of oscillating test functions, Γ-convergence, two-scale convergence, and the unfolding method. We refer to [12] for an accessible introduction to the mathematics of homogenization and for an overview of the different homogenization methods.
We are aware of two earlier works [4,34] containing rigorous homogenization results for the bidomain model (but see [15,16,43] for examples of elliptic and parabolic equations on"two-component" domains). With a fairly advanced proof involving Γ-convergence, the De Giorgi "minimizing movement" approach, timediscretization, variational problems, and two limit procedures, the homogenization result in [34] covers the generalized FitzHugh-Nagumo ionic model [17]. The proof of the result in [4] is more basic in the sense that it employs only two-scale convergence arguments, but it handles only a restricted class of ionic models.
We mention that there are several complications preventing the application of standard homogenization results (for elliptic/parabolic equations) to the bidomain equations, including its degenerate structure (seen at the tissue-level), resulting from differing anisotropies of the intra-and extracellular spaces, and the highly nonlinear, oscillating dynamic boundary condition (seen at the cell-level).
The main contribution of our paper is to provide a simple homogenization proof that can handle some relevant nonlinear membrane models (the generalized FitzHugh-Nagumo model), relying only on basic two-scale convergence techniques.
We now explain our contribution in more detail. The point of departure is the following microscopic model [13,14,22,41] for the electric activity in cardiac tissue: ε (∂ t v ε + I(v ε , w ε )) = −ν · σ ε e ∇u ε e on (0, T ) × Γ ε , ∂ t w ε = H(v ε , w ε ) on (0, T ) × Γ ε , (1.1) where ν denotes the unit normal pointing out of Ω ε i (and into Ω ε e ). Cardiac tissue consists of an assembly of elongated cylindrical-shaped cells coupled together (endto-end and side-to-side) to provide intercellular communication. The entire cardiac domain Ω ⊂ R 3 is viewed as a "two-component" domain and split into two εperiodic open sets Ω ε i , Ω ε e corresponding to the intra-and extracellular spaces. The sets Ω ε i , Ω ε e , which are assumed to be disjoint and connected, are separated by an ε-periodic surface Γ ε representing the cell membrane, so that Ω = Ω ε i ∪ Ω ε e ∪ Γ ε . The main geometrical assumption is that the intra-and extracellular domains are ε-dilations of some reference cells Y i , Y e ⊂ Y := [0, 1] 3 , periodically repeated over R 3 . Although our results are valid for general Lipschitz domains, for simplicity of presentation, we assume that the entire cardiac domain Ω is the open cube In (1.1), σ ε j is the conductivity tensor and s ε j is the stimulation current, relative to Ω ε j for j = i, e. The functions s ε i , s ε e are assumed to be at least bounded in L 2 , independently of ε. As usual in homogenization theory, the conductivity tensors σ ε i , σ ε e are assumed to have the form where σ j = σ j (x, y) satisfies the usual conditions of uniform ellipticity and periodicity (in y). Despite the fact that the inhomogeneities of the domains impose ε-oscillations in the conductivity tensors (via gap junctions), the main source of inhomogeneity in the microscopic model is not the conductivities σ ε i and σ ε e , but the domains Ω ε i and Ω ε e themselves. We allow for inhomogeneous and oscillating conductivities for the sake of generality.
We denote by u ε j the electric potential in Ω ε j (j = i, e). On Γ ε , v ε := u ε i − u ε e is the transmembrane potential and I(v ε , w ε ) is the ionic current depending on v ε and a gating variable w ε . The left-hand side of the third and fourth equations in (1.1) describes the current across the membrane as having a capacitive component, depending on the time derivative of the transmembrane potential, and a nonlinear ionic component I corresponding to the chosen membrane model. In this article we consider the generalized FitzHugh-Nagumo model [17]. We choose to focus on this membrane model for definiteness, but our arguments can be adapted to many other models satisfying reasonable technical assumptions [7,8,13,38,41,42]. For each fixed ε > 0, the functions σ ε j , s ε j , I, H in (1.1) are given and we wish to solve for (u ε i , u ε e , v ε , w ε ). To this end, we must augment the system (1.1) with initial conditions for v ε , w ε and Neumann-type boundary conditions for u ε i , u ε e (ensuring no current flow out of the heart): where n is the outward unit normal to Ω. It is proved in [14,41] that the microscopic bidomain model (1.1), (1.3) possesses a unique weak solution. This solution satisfies a series of a priori estimates. For us it is essential to know how these estimates depend on the parameter ε. We will therefore outline a proof of these estimates. The dimensionless number ε is a small positive number representing the ratio of the microscopic and macroscopic scales, that is, considering Ω as fixed, it is proportional to the cell diameter. The goal of homogenization is to investigate the limit of a sequence of solutions {(u ε i , u ε e , v ε , w ε )} ε>0 to (1.1), (1.3). By the multiple scales method [12,14,22], the electric potentials u ε i , u ε e , v ε and the state variable w ε exhibit the following asymptotic expansions in powers of the parameter ε: where y = x/ε denotes the microscopic variable, and each term in the expansions is a function of both the slow (macroscopic) variable x and the fast (microscopic) variable y, periodic in y. Substituting the above expansions into (1.1), and equating all terms of the same orders in powers of ε, we obtain after some routine arguments that the zero order terms u i , u e , v, w are independent of the fast variable y and satisfy the (macroscopic) bidomain model [14,22,34] (1.4) where the homogenized conductivity tensors M i (x), M e (x) are given by and the y-periodic (vector-valued) function χ j = χ j (x, y) solves the cell problem on Ω × Γ.
Note that the effective potentials u i , u e in (1.4) are defined at every point of Ω, while in the microscopic model they live on disjoint sets Ω ε i , Ω ε e . In (1.4), (1.5), (1.6) the sets Y i , Y e are the intra and extracellular spaces within the reference unit cell Y , separated by the cell membrane Γ (see Section 2 for details). It is worth noting that the bidomain model is often stated in terms of the "geometric" parameter χ = |Γ| |Y | = |Γ| representing the surface-to-volume ratio of the cardiac cells. Regarding the existence and uniqueness of properly defined solutions, standard theory for parabolic-elliptic systems does not apply naturally to the bidomain model (1.4). A number of works [6,7,8,14,42] have recently provided well-posedness results for (1.4), applying differing solution concepts and technical frameworks.
As alluded to earlier, we will provide a rigorous derivation of the homogenized system (1.4), (1.5), (1.6) based on the theory of two-scale convergence (see [33] and [1,2,28]). This result is not covered by standard parabolic homogenization theory. A complication is the nonlinear dynamic boundary condition (posed on an underlying oscillating surface), which makes it difficult to pass to the limit in (1.1) as ε → 0. The aim is to prove that a sequence {u ε i , u ε e , v ε , w ε } ε>0 of solutions to the microscopic problem two-scale converges to the solution of the bidomain model (1.4). However, two-scale convergence is not "strong enough" to justify passing to the limit in the nonlinear boundary condition. To handle this difficulty we use the boundary unfolding operator [11], establishing strong convergence of T b ε (v ε ) in L 2 ((0, T ) × Ω × Γ), where T b ε denotes the boundary unfolding operator. The boundary unfolding operator makes our proof flexible enough to handle a range of membrane models, exemplified by the generalized FitzHugh-Nagumo model. Unfolding operators, presented and carefully analyzed in [10,11], facilitate elementary proofs of classical homogenization results on fixed as well as perforated domains/surfaces. An unfolding operator T ε maps a function v(x) defined on an oscillating domain/surface to a higher dimensional function T ε (v)(x, y) on a fixed domain, to which one can apply standard convergence theorems in fixed L p spaces. Reflecting the "two-component" nature of the cardiac domain, it makes sense to use two unfolding operators T i,e ε , linked to the intra-and extracellular domains Ω ε i,e . In this paper, however, we mainly unfold functions defined on the cell membrane, utilizing the boundary unfolding operator T b ε . For somewhat similar unfolding of "two-component" domains separated by a periodic boundary, see [15,16,43]. For other relevant works that combine twoscale convergence and unfolding methods, we refer to [19,18,20,29,32]. Among these, our work borrows ideas mostly from [19,18,32].
The remaining part of the paper is organized as follows: In Section 2, we collect relevant functional spaces and analysis results. Moreover, we gather definitions and tools linked to two-scale convergence and unfolding operators. In Section 3, we define precisely what is meant by a weak solution of the microscopic problem (1.1), state a well-posedness result, and establish several "ε-independent" a priori estimates. The main homogenization result is stated and proved in Section 4.

Preliminaries
2.1. Some functional spaces and tools. For a general review of integer and fractional order Sobolev spaces (on Lipschitz domains) and relevant analysis tools, see [9,Chaps. 2 & 3] and [30,Chap. 3]. For relevant background material on mathematical homogenization, we refer to [12].
Let Ω ⊂ R 3 be a bounded open set with Lipschitz boundary. We denote by C ∞ 0 (Ω) the (infinitely) smooth functions with compact support in Ω. The space of smooth Y -periodic functions is denoted by C ∞ per (Y ). The closure of this space under the norm ∇(·) L 2 (Y ) is denoted by H 1 per (Y ). We write H s for the L 2 -based Sobolev spaces W s,2 (s ∈ (0, 1]).
We make use of Sobolev spaces on surfaces, as defined for example in [27, p. 34] and [30, p. 96]. Specifically, we use the (Hilbert) space H 1/2 (Γ), for a two-dimensional Lipschitz surface Γ ⊂ Ω, equipped with the norm and dS is the two-dimensional surface measure. We define the dual space of H 1/2 (Γ) as H −1/2 (Γ) := (H 1/2 (Γ)) * , equipped with the norm of dual spaces The following trace inequality holds: Any function in H 1/2 (Γ) can be characterized as the trace of a function in H 1 (Ω).
The trace map has a continuous right inverse I : where the constant C depends only on Γ. We need the Sobolev inequality , 4]. This embedding is compact for p ∈ [1, 4). In particular, H 1/2 (Γ) is compactly embedded in L 2 (Γ). Let X be a separable Banach space X and p ∈ [1, ∞], We make routinely use of Lebesgue-Bochner spaces such as L p (Ω; X) and L p (0, T ; X). We also use the spaces of continuous functions from Ω to X and (0, T ) to X, denoted by C(Ω; X) and C(0, T ; X), respectively, and the similar spaces with C replaced by C p or C p 0 . If X is a Banach space, then X/R denotes the space consisting of classes of functions in X that are equal up to an additive constant.
Let us dwell a bit further on the time continuity of functions in V Γ,T . By (2.3), H 1/2 ⊂ L 4 (Γ) and so L 4/3 (Γ) ⊂ H −1/2 (Γ). Therefore, Integrating this inequality over t ∈ (0, T − ∆ t ), accompanied by a few elementary manipulations, results in the temporal translation estimate There is a compact embedding of V Γ,T in L 2 (0, T ; L 2 (Γ)). As pointed out above, We need a generalization of this result due to Simon [37]. Given two Banach spaces X 1 ⊂ X 0 , with X 1 compactly embedded in X 0 , let K be a collection of functions in L p (0, T ; X 0 ), p ∈ [1, ∞]. The work [37] supplies several results ensuring the compactness of K in L p (0, T ; X 0 ) (in C([0, T ]; X 0 ) if p = ∞). For example, we can assume that K is bounded in L 1 loc (0, T ; X 1 ) and u(· + ∆ t ) − u L p (0,T −∆t;X0) → 0 as ∆ t → 0, uniformly for u ∈ K, cf. [37,Theorem 3]. We apply this result with p = 2, Denote by K ε the set of k ∈ Z 3 for which k∈K ε ε (k + Y ) = Ω. We define the intracellular domain Ω ε i , the extracellular domain Ω ε e , and the cell membrane Γ ε as Both sets Ω ε i , Ω ε e are connected Lipschitz domains, see Figure 1.1. Note however that it is impossible to have both Ω ε i and Ω ε e connected in a two-dimensional picture. Figure 1. The rescaled sets Ω ε i , Ω ε e , Γ ε (left) and the unit cell Y (right).
To derive estimates for the microscopic model, we employ the following trace inequality for ε-periodic hypersurfaces: We need a uniform Poincaré inequality for perforated domains [11].
for all u ∈ H 1 Ω ε j , j = i, e. Estimate (2.9) holds under mild regularity assumptions on the perforated domains; a Lipschitz boundary is more than sufficient (but connectedness is essential).
Recall that a sequence . We express this symbolically as By density properties, the convergence (2.10) also holds for test functions ϕ from The two-scale compactness theorem [1,12] is of fundamental importance.
Consider a sequence {u ε j } ε>0 of functions defined on the perforated domain (0, T ) × Ω ε j , j = i, e. We write u ε j for the zero-extension of u ε j to (0, T ) × Ω: has a two-scale convergent subsequence, provided we know that u ε j L 2 ((0,T )×Ω ε j ) ≤ C ∀ε > 0. However, this is not true in general for the gradient of u ε j , even if u ε j L 2 (0,T ;H 1 (Ω ε j )) ≤ C, since the extension by zero creates a discontinuity across Γ ε . Instead the following statement holds true for the gradient: Then there exist a subsequence ε n → 0 and functions u ∈ L 2 (0, T ; H 1 (Ω)), Here, ½ Yj (y) denotes the characteristic function of Y j , For a proof of this lemma in the time independent case, see [1, Theorem 2.9]. The extension to time dependent functions is straightforward.
There is an extension of two-scale convergence to periodic surfaces [2]. Recall that a periodic surface Γ ε is given by As with (2.10), this convergence continues to hold for test functions ϕ that belong to L 2 ((0, T ) × Ω; C per (Γ)).
There is a version of Theorem 2.2 for functions on periodic surfaces [2].
for some function C that is independent of ε > 0. Then there exist a subsequence ε n → 0 and a function v ∈ L 2 ((0, T ) × Ω; L 2 (Γ)) such that as n → ∞, One can characterize the two-scale limit of traces of bounded sequences in L 2 (0, T ; H 1 (Ω ε )) as the trace of the two-scale limit [2].
Applying the above decomposition to x/ε gives The advantage of the unfolding operator is that we can formulate questions of convergence in a fixed space L 2 (0, T ; L 2 (Ω × Γ)). All definitions and results in this section are formulated in L 2 spaces. Everything remains the same, however, if we replace L 2 by L p for any p ∈ [1, ∞). We refer to [11] for the definition of the boundary unfolding operator and proofs of the properties listed next.

Microscopic bidomain model
In this section we present a relevant notion of (weak) solution for the microscopic problem (1.1), (1.3), along with an accompanying existence theorem. We also derive some "ε-independent" a priori estimates, which are used later to extract two-scale convergent subsequences.
3.1. Assumptions on the data. We impose the following set of assumptions on the "membrane" functions I, H: • Generalized FitzHugh-Nagumo model: For v, w ∈ R, for some constants c I > 0 and c I,1 , c I,2 , c I,3 , c I,4 , c H,2 ≥ 0. The classical FitzHugh-Nagumo model corresponds to where a ∈ (0, 1) and k, ǫ > 0 are constants. Repeated applications of Cauchy's inequality yields for some constants γ > 0 and β ≥ 0. This inequality will be used to bound the transmembrane potential in the L 4 norm. Consider a quadratic matrix A, which always can be written as the sum of its symmetric part 1 2 (A+A ⊤ ) and its skew-symmetric part 1 2 (A−A ⊤ ). Recall that in a quadratic form z → Az · z the skew-symmetric part does not contribute. Therefore, letting λ min and λ max denote respectively the minimum and maximum eigenvalues of the symmetric part of A, we have λ min |z| 2 ≤ Az · z ≤ λ max |z| 2 , ∀z.
For µ > 0, consider the function F µ : R 2 → R 2 defined by Denote by λ µ min (z), λ µ max (z) the minimum, maximum eigenvalues of the symmetric part of the matrix ∇F µ (z). To ensure that weak solutions are unique, we need an additional assumption on I, H expressed via F µ [8, p. 479]: ∃µ, λ > 0 such that

One can verify that the FitzHugh-Nagumo model (3.1) obeys (3.3) (with µ = εk).
A consequence of (3.3) is that ∇F µ (z)z · z ≥ λ |z| 2 for allz, z ∈ R 2 . Therefore, writing More explicitly, assumption (3.3) implies the following "dissipative structure" on a suitable linear combination of I and H: This inequality implies the L 2 stability (and thus uniqueness) of weak solutions.
Remark 3.1. There are many membrane models of cardiac cells [13,38]. We utilize the FitzHugh-Nagumo model [17], which is a simplification of the Hodgin-Huxley model of voltage-gated ion channels. It is possible to treat other membrane models by blending the arguments used herein with those found in [7,8,13,38,41,42].
As a natural assumption for homogenization, we assume that the ε-dependence of the conductivities σ ε j (j = i, e), the applied currents s ε j (j = i, e), and the initial data v ε 0 , w ε 0 decouples into a "fast" and a "slow" variable: for some fixed functions σ j (x, y), s j (x, y), v 0 (x, y), w 0 (x, y) that are Y -periodic in the second argument. The conductivity tensors are assumed to be bounded and continuous, and satisfy the usual ellipticity condition, i.e., there exists α > 0 such that for j = i, e. Finally, we assume that each σ j is symmetric: σ T j = σ j . The regularity assumption (3.5) implies that σ j is an admissible test function for two-scale convergence [1], which means that for every two-scale convergent sequence {ϕ ε } ε>0 , ϕ ε 2 ⇀ ϕ. This convergence still holds if the second part of (3.5) is replaced by σ j ∈ L 2 (Ω; C per (Y j )).
For the stimulation currents we assume the compatibility condition j=i,e Ω ε j s ε j dx = 0, and the boundedness in L 2 : which, in view of (3.4), is guaranteed if we take [12, p. 174] (3.9) s j ∈ L 2 (Ω; C per (Y j )).
Throughout this paper we denote by C a generic constant, not depending on the parameter ε. The actual value of C may change from one line to the next.

Weak solutions.
Testing (1.1) against appropriate functions ϕ i , ϕ e , ϕ w we obtain the weak formulation of the microscopic bidomain model (1.1), (1.3), cf. [14,41] for details. We note that the terms involving the boundary ∂Ω vanish due to the Neumann boundary condition (1.3).
Definition 3.2 (weak formulation of microscopic system). A weak solution to (1.1), (1.3) is a collection (u ε i , u ε e , v ε , w ε ) of functions satisfying the following conditions: i (algebraic relation).
Remark 3.4. Consider a weak solution (u ε i , u ε e , v ε , w ε ) according to Definition 3.2. Thanks to (2.1), the trace of ϕ j ∈ L 2 (0, T ; H 1 (Ω ε j )) belongs to L 2 (0, T ; H 1/2 (Γ ε )). By the Sobolev inequality (2.3), the trace of ϕ j belongs also to L 2 (0, T ; L 4 (Γ ε )). In Definition 3.2 we ask additionally that the trace of ϕ j belongs to L 4 ((0, T ) × Γ ε ) to ensure that the surface terms in (3.13), (3.14) are well-defined Indeed, T 0 ∂ t v ε , ϕ j dt is well-defined for such ϕ j . Moreover, For the membrane model (GFHN), the growth condition on I implies and therefore I(v ε , w ε ) Consequently, J < ∞. We actually have a more precise bound. As Integrating this over t ∈ (0, T ) yields The integral on the right-hand side of (3.15) can be treated similarly, since The remaining integrals are trivially well defined.

3.3.
Existence of solution and a priori estimates. Existence and uniqueness results for certain classes of membrane models have been established in [14,41]. These works employ the variable v = u i − u e to convert (1.1) into a non-degenerate "abstract" parabolic equation. The authors in [14] then appeal to the theory of variational inequalities, whereas in [41] the Schauder fixed point theorem is applied to conclude the existence of a solution.
The following theorem can be proved by adapting arguments found in [14,41] (see also [21]), or those utilized in [8,6,5] for the (macroscopic) bidomain model.  The next lemma, which is utilized below to derive some a priori estimates, is a consequence of the uniform Poincaré inequality (2.9) and the trace inequality for ε-periodic surfaces (2.8). A similar result is used in [34].
Lemma 3.7. Let u j ∈ H 1 (Ω ε j ), j = i, e, Ω ε e u e dx = 0, and set v := u i Γ ε − u e Γ ε . There is a positive constants C, independent of ε, such that Proof. First, since u e = 0, To estimate the L 2 norm of u i , write u i =ū i +ũ i , whereū i : In view of the Poincaré inequality (2.9), we obtain , where the final inequality is a result of (3.18) and (3.19).
For the sake of the upcoming homogenization result, we now list some precise (ε-independent) a priori estimates.
Proof. We only outline a proof of these (mostly standard) estimates.
For the upcoming convergence analysis, we need a temporal translation estimate for the membrane potential v ε . Lemma 3.10 (temporal translation estimate in L 2 ). Referring back to Theorem 3.5, let (u ε i , u ε e , v ε , w ε ) be the weak solution of (1.1), (1.3). Then v ε satisfies the following L 2 temporal translation estimate for some ε-independent constant C > 0, for any sufficiently small ∆ t > 0: Proof. The translated functions (u ε i , u ε e , v ε , w ε )(t + ∆ t ) constitute a weak solution of (1.1) on (0, T − ∆ t ) with initial data (v ε , w ε )(∆ t ) and stimulation currents (s ε i , s ε e )(t + ∆ t ). We subtract the equations (linked to Ω ε i and Ω ε e ) for the original weak solution from the equations satisfied by the translated one, and add the resulting equations. The result is Specifying the test functions ϕ i , ϕ e , ϕ w as Let us write this equation as By the Cauchy-Schwarz and Minkowski integral inequalities, Similarly, |J 2 | ≤ C 2 ∆ t . We need the following facts involving a real-valued function f ∈ L p (p > 1): where the second line follows from Young's convolution inequality. By (3.5) and the Cauchy-Schwarz inequality, .

As a result of Minkowski integral inequality,
An application of Hölder's inequality yields where we have repeated the argument for J 1 involving (3.26), with the implication that . Moreover, we have used the basic inequality ab ≤ 3 4 a 4/3 + 1 4 b 4 for positive numbers a, b. Similarly, and Summarizing our findings, we conclude that (3.25) holds.
Remark 3.11. Due to the degenerate structure of the microscopic system (1.1), temporal estimates are not available for the intra-and extracellular potentials u ε i , u ε e .

The homogenization result
This section contains the main result of the paper. We start by recalling the weak formulation of the macroscopic bidomain system (1.4), which is augmented with the following initial and boundary conditions: n · σ∇u j = 0 on (0, T ) × ∂Ω, j = i, e. ii (regularity).
The macroscopic bidomain system is well studied for a variety of cellular models [5,6,7,8,14,42]. For the following result, see [5,6,7]. We are now in a position to state the main result, which should be compared to Theorem 1.3 in [34].
Remark 4.4. The strong convergence (4.3) is a corrector-type result. In the current setting it allows us to pass to the limit in the nonlinear ionic terms. By employing standard techniques [12] one can also show that the energies converge to the homogenized energy The rest of this section is devoted to the proof of Theorem 4.3. Homogenization of the linear terms in (1.1) is handled with standard techniques, cf. Subsection 4.1.
Passing to the limit in the nonlinear terms I, H is more challenging. Although the intracellular/extracellular functions u ε i (defined on Ω ε i ) and u ε e (defined on Ω ε e ) do not converge strongly, some kind of strong compactness is expected for the εscaled version of the transmembrane potential v ε = u ε i Γ ε − u ε e Γ ε (defined on Γ ε ), since we control both the temporal (3.25) and spatial (fractional) derivatives (3.20).
Wild oscillations of the underlying domain do however pose difficulties. For this reason, we use the boundary unfolding operator T b ε , cf. (2.12), to transform the problem of convergence on the oscillating set Γ ε to the fixed set Ω × Γ. Roughly speaking, (3.20) is used to conclude that T b ε (v ε ) is uniformly bounded in L 2 t L 2 x H 1/2 y .
Nonlinear functions are not continuous with respect to weak convergence, which prevents us from immediately making the identifications I = I(v, w), H = H(v, w).
Using the two-scale convergences in (4.4), Remark 3.3, the choice of test function ϕ(t, x) + εϕ 1 t, x, x ε , with ϕ ∈ C ∞ 0 ((0, T ) × Ω), ϕ 1 ∈ C ∞ 0 ((0, T ) × Ω; C ∞ per (Y )), in (3.13), (3.14), and (3.15), standard manipulations [12,2] will reveal that the two-scale limit (u i , u e , v, w) satisfies the following equations:  where v, w are identified in (4.4). Besides, since w ε appears linearly in I and H, cf. (GFHN), the weak convergence of T b ε (w ε ) ε>0 in L 2 ((0, T ) × Ω × Γ) is enough to pass to the limit in (4.8), if we establish strong convergence of T b ε (v ε ) ε>0 . As a step towards verifying the required strong convergence, we need to show that (t, y) → T b ε (v ε ) ε>0 is strongly precompact in L 2 t,y , for fixed x ∈ Ω. As a result of Lemma 3.8, T b ε (v) is bounded in L 2 t L 2 x H 1/2 y , uniformly in ε. However, according to Lemma 3.8, the time derivative ∂ t v ε is merely of order 1/ε in the L 2 t H −1/2 x norm. Therefore, we cannot expect ∂ t T b ε (v ε ) (assuming that this object is meaningful) to be bounded in L 2 t L 2 x H −1/2 y , uniformly in ε. As a consequence, strong L 2 t,y compactness of T b ε (v ε ) ε>0 is not deducible from the classical Aubin-Lions theorem. Instead of attempting to control (in a negative space) the whole derivative ∂ t T b ε (v ε ), we will make use of a temporal translation estimate with respect to the L 2 norm (cf. lemma below). The L 2 t,y compactness will then be a consequence of the Simon lemma (cf. Subsection 2.1).
The rest of this section is devoted to the detailed proof that T b ε (v ε ) ε>0 is strongly precompact in the fixed space L 2 ((0, T ) × Ω × Γ). We begin with where we have also used εY k,m because w ε 2−S ⇀ w, cf. (4.4) and (2.17). Besides, w ε appears linearly in I and H, cf. (GFHN). We have shown that T b ε (v ε ) ε>0 is strongly precompact in L 2 ((0, T ) × Ω × Γ). It then follows that T b ε (v ε ) → v in L 2 ((0, T ) × Ω × Γ) and a.e. in (0, T ) × Ω × Γ, along a subsequence as ε → 0 (not relabelled), where v is the two-scale limit of {v ε } ε>0 identified in (4.4). By way of estimate (d) in (3.20) and (the L p version of) (2.15), where C is independent of ε. In view of this estimate and the Vitali convergence theorem, we conclude the validity of (4.8). This finishes the proof of Theorem 4.3.