The cardiac bidomain model and homogenization

1.   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 $ \sigma_i,\sigma_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 $ \sigma_i,\sigma_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, $ \Gamma $-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 [16,15,43] for examples of elliptic and parabolic equations on"two-component" domains). With a fairly advanced proof involving $ \Gamma $-convergence, the De Giorgi "minimizing movement" approach, time-discretization, 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:

where $ \nu $ denotes the unit normal pointing out of $ \Omega_i^{\varepsilon} $ (and into $ \Omega_e^{\varepsilon} $). Cardiac tissue consists of an assembly of elongated cylindrical-shaped cells coupled together (end-to-end and side-to-side) to provide intercellular communication. The entire cardiac domain $ \Omega\subset \mathbb{R} ^3 $ is viewed as a "two-component" domain and split into two $ \varepsilon $-periodic open sets $ \Omega_i^{\varepsilon} $, $ \Omega_e^{\varepsilon} $ corresponding to the intra- and extracellular spaces. The sets $ \Omega_i^{\varepsilon} , \Omega_e^{\varepsilon} $, which are assumed to be disjoint and connected, are separated by an $ \varepsilon $-periodic surface $ \Gamma^{\varepsilon} $ representing the cell membrane, so that $ \Omega = \Omega_i^{\varepsilon} \cup \Omega_e^{\varepsilon} \cup \Gamma^{\varepsilon} $. The main geometrical assumption is that the intra- and extracellular domains are $ \varepsilon $-dilations of some reference cells $ Y_i,Y_e\subset Y: = [0,1]^3 $, periodically repeated over $ \mathbb{R} ^3 $. Although our results are valid for general Lipschitz domains, for simplicity of presentation, we assume that the entire cardiac domain $ \Omega $ is the open cube

In (1.1), $ \sigma_j^{ \varepsilon } $ is the conductivity tensor and $ s_j^{ \varepsilon } $ is the stimulation current, relative to $ \Omega_j^{\varepsilon} $ for $ j = i,e $. The functions $ s_i^ \varepsilon ,s_e^ \varepsilon $ are assumed to be at least bounded in $ L^2 $, independently of $ \varepsilon $. As usual in homogenization theory, the conductivity tensors $ \sigma_i^{ \varepsilon }, \sigma_e^{ \varepsilon } $ are assumed to have the form

where $ \sigma_j = \sigma_j(x,y) $ satisfies the usual conditions of uniform ellipticity and periodicity (in $ y $). Despite the fact that the inhomogeneities of the domains impose $ \varepsilon $-oscillations in the conductivity tensors (via gap junctions), the main source of inhomogeneity in the microscopic model is not the conductivities $ \sigma_i^ \varepsilon $ and $ \sigma_e^ \varepsilon $, but the domains $ \Omega_i^{\varepsilon} $ and $ \Omega_e^{\varepsilon} $ themselves. We allow for inhomogeneous and oscillating conductivities for the sake of generality.

We denote by $ u_j^{ \varepsilon } $ the electric potential in $ \Omega_j^{\varepsilon} $ ($ j = i,e $). On $ \Gamma^{\varepsilon} $, $ v^{ \varepsilon } : = u_i^{ \varepsilon }-u_e^{ \varepsilon } $ is the transmembrane potential and $ I(v^{ \varepsilon },w^{ \varepsilon }) $ is the ionic current depending on $ v^{ \varepsilon } $ and a gating variable $ w^{ \varepsilon } $. 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 $ \varepsilon >0 $, the functions $ \sigma_j^{ \varepsilon }, s_j^{ \varepsilon }, I,H $ in (1.1) are given and we wish to solve for $ (u_i^ \varepsilon ,u_e^ \varepsilon ,v^ \varepsilon ,w^ \varepsilon ) $. To this end, we must augment the system (1.1) with initial conditions for $ v^ \varepsilon , w^ \varepsilon $ and Neumann-type boundary conditions for $ u_i^ \varepsilon , u_e^ \varepsilon $ (ensuring no current flow out of the heart):

where $ n $ is the outward unit normal to $ \Omega $. 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 $ \varepsilon $. We will therefore outline a proof of these estimates.

The dimensionless number $ \varepsilon $ is a small positive number representing the ratio of the microscopic and macroscopic scales, that is, considering $ \Omega $ as fixed, it is proportional to the cell diameter. The goal of homogenization is to investigate the limit of a sequence of solutions $ \left\{\left(u^{ \varepsilon }_i,u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon }\right)\right\}_{ \varepsilon >0} $ to (1.1), (1.3). By the multiple scales method [12,14,22], the electric potentials $ u^{ \varepsilon }_i,u_e^{ \varepsilon },v^{ \varepsilon } $ and the state variable $ w^ \varepsilon $ exhibit the following asymptotic expansions in powers of the parameter $ \varepsilon $:

where $ y = x/ \varepsilon $ 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 $ \varepsilon $, 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]

where the homogenized conductivity tensors $ M_i(x), M_e(x) $ are given by

and the $ y $-periodic (vector-valued) function $ \chi_j = \chi_j(x,y) $ solves the cell problem

Note that the effective potentials $ u_i, u_e $ in (1.4) are defined at every point of $ \Omega $, while in the microscopic model they live on disjoint sets $ \Omega_i^{\varepsilon} , \Omega_e^{\varepsilon} $. 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 $ \Gamma $ (see Section 2 for details). It is worth noting that the bidomain model is often stated in terms of the "geometric" parameter $ \chi = \frac{\left|\Gamma\right|}{\left|Y\right|} = \left|\Gamma\right| $ 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 $ \varepsilon \to 0 $. The aim is to prove that a sequence $ \left\{u_i^{ \varepsilon },u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon }\right\}_{ \varepsilon >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 [10], establishing strong convergence of $ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon }) $ in $ L^2((0,T)\times \Omega \times \Gamma ) $, where $ \mathcal{T} _{ \varepsilon }^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 [11,10], facilitate elementary proofs of classical homogenization results on fixed as well as perforated domains/surfaces. An unfolding operator $ \mathcal{T} _ \varepsilon $ maps a function $ v(x) $ defined on an oscillating domain/surface to a higher dimensional function $ \mathcal{T} _{ \varepsilon }(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 $ \mathcal{T} _ \varepsilon ^{i,e} $, linked to the intra- and extracellular domains $ \Omega ^ \varepsilon _{i,e} $. In this paper, however, we mainly unfold functions defined on the cell membrane, utilizing the boundary unfolding operator $ \mathcal{T} _ \varepsilon ^b $.

For somewhat similar unfolding of "two-component" domains separated by a periodic boundary, see [16,15,43]. For other relevant works that combine two-scale 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 "$ \varepsilon $-independent" a priori estimates. The main homogenization result is stated and proved in Section 4.

2.   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 $ \Omega \subset \mathbb{R} ^3 $ be a bounded open set with Lipschitz boundary. We denote by $ C_0^{\infty}( \Omega ) $ the (infinitely) smooth functions with compact support in $ \Omega $. The space of smooth $ Y $-periodic functions is denoted by $ C^{\infty}_{{\rm{per}}}(Y) $. The closure of this space under the norm $ \|\nabla(\cdot)\|_{L^2(Y)} $ is denoted by $ H_{{\rm{per}}}^1(Y) $. We write $ H^s $ for the $ L^2 $-based Sobolev spaces $ W^{s,2} $ ($ s\in (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}( \Gamma ) $, for a two-dimensional Lipschitz surface $ \Gamma \subset \Omega $, equipped with the norm

where

and $ dS $ is the two-dimensional surface measure. We define the dual space of $ H^{1/2}( \Gamma ) $ as $ H^{-1/2}( \Gamma ): = (H^{1/2}( \Gamma ))^* $, equipped with the norm of dual spaces

The following trace inequality holds:

Any function in $ H^{1/2}( \Gamma ) $ can be characterized as the trace of a function in $ H^1( \Omega ) $. The trace map has a continuous right inverse $ \mathcal{I}:H^{1/2}( \Gamma )\rightarrow H^1( \Omega ) $, satisfying

where the constant $ C $ depends only on $ \Gamma $. We need the Sobolev inequality

Indeed, $ H^{1/2}( \Gamma ) $ is continuously embedded in $ L^p( \Gamma ) $ for $ p\in [1,4] $. This embedding is compact for $ p\in [1,4) $. In particular, $ H^{1/2}( \Gamma ) $ is compactly embedded in $ L^2( \Gamma ) $.

Let $ X $ be a separable Banach space $ X $ and $ p\in [1,\infty] $, We make routinely use of Lebesgue-Bochner spaces such as $ L^p(\Omega;X) $ and $ L^p(0,T;X) $. We also use the spaces of continuous functions from $ \Omega $ to $ X $ and $ (0,T) $ to $ X $, denoted by $ C( \Omega ;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/ \mathbb{R} $ denotes the space consisting of classes of functions in $ X $ that are equal up to an additive constant.

Recall that $ H^{1/2}( \Gamma ) $ is a Hilbert space embedded in a continuous and dense way in $ L^2( \Gamma ) $. The (Lions-Magenes) integration-by-parts formula holds for functions $ u_1,u_2 $ that belong to the Banach space

equipped with the norm

where $ \left\|u\right\|_{X_1\cap X_2} = \max \left(\left\|u\right\|_{X_1},\left\|u\right\|_{X_2} \right) $, $ \left\|u\right\|_{X_1+X_2} = \inf\limits_{u = u_1+u_2} \left(\left\|u_1\right\|_{X_1} +\left\|u_2\right\|_{X_2}\right) $. Indeed, for $ u_1,u_2 \in \mathcal{V}_{ \Gamma ,T} $, $ [0,T]\ni t\mapsto \left ( u_1(t), u_2(t) \right)_{L^2( \Gamma )} $ is continuous and

for all $ t_1,t_2\in [0,T] $, $ t_1<t_2 $, where $ \left\langle \cdot, \cdot \right \rangle $ denotes the duality pairing between $ H^{-1/2}( \Gamma )+L^{4/3}( \Gamma ) $ and $ H^{1/2}( \Gamma )\cap L^4( \Gamma ) $. For a proof of (2.4) that can be adapted to our situation, see e.g. [9,p. 99].

Taking $ u_1 = u_2 = u\in \mathcal{V}_{ \Gamma ,T} $, we obtain the chain rule

Adapting standard arguments (see e.g. [9,p. 101]), the embedding

is continuous. Indeed, this result follows from the continuity of the squared norm $ t\mapsto \left\|u(t)\right\|_{L^2( \Gamma )}^2 $ (see above) and the weak continuity of $ u $ in $ L^2( \Gamma ) $. The latter results from an easily obtained bound on $ u $ in $ L^\infty(0,T;L^2( \Gamma )) $ and the continuity of $ u $ in $ H^{-1/2}( \Gamma ) $, both facts being deducible from (2.4).

Let us dwell a bit further on the time continuity of functions in $ \mathcal{V}_{ \Gamma ,T} $. By (2.3), $ H^{1/2}\subset L^4( \Gamma ) $ and so $ L^{4/3}( \Gamma )\subset H^{-1/2}( \Gamma ) $. Therefore,

With $ u_1 = u\in \mathcal{V}_{ \Gamma ,T} $ and $ u_2 = \phi\in H^{1/2}( \Gamma ) $ in (2.4), it follows that

Fix a small shift $ \Delta_t>0 $. Specifying $ t_1 = t\in (0,T-\Delta_t) $, $ t_2 = t+\Delta_t $, and $ \phi = u(t+\Delta_t,\cdot)-u(t,\cdot) $ gives

Integrating this inequality over $ t\in (0,T-\Delta_t) $, accompanied by a few elementary manipulations, results in the temporal translation estimate

where $ C_T = 2^{1/2}T^{1/4} $. A similar estimate holds for negative $ \Delta_t $.

There is a compact embedding of $ \mathcal{V}_{ \Gamma ,T} $ in $ L^2(0,T;L^2( \Gamma )) $. As pointed out above, $ \mathcal{V}_{ \Gamma ,T} $ is a subset of $ \left\{u\in L^2(0,T;H^{1/2}( \Gamma ): \,{\partial _t}u \in L^{4/3}(0,T;H^{-1/2}( \Gamma ))\right\} $, which is compactly embedded in $ L^2(0,T;L^2( \Gamma )) $ by the Aubin-Lions theorem.

We need a generalization of this result due to Simon [37]. Given two Banach spaces $ X_1\subset X_0 $, with $ X_1 $ compactly embedded in $ X_0 $, let $ \mathcal{K} $ be a collection of functions in $ L^p(0,T;X_0) $, $ p\in [1,\infty] $. The work [37] supplies several results ensuring the compactness of $ \mathcal{K} $ in $ L^p(0,T;X_0) $ (in $ C([0,T];X_0) $ if $ p = \infty $). For example, we can assume that $ \mathcal{K} $ is bounded in $ L^1_{{\rm{loc}}}(0,T;X_1) $ and

cf. [37,Theorem 3]. We apply this result with $ p = 2 $, $ X_1 = H^{1/2}( \Gamma^{\varepsilon} ) $, $ X_0 = L^2( \Gamma^{\varepsilon} ) $. Another result involves a third Banach space $ X_{-1} $ (e.g. $ X_{-1} = H^{-1/2}( \Gamma^{\varepsilon} ) $), such that $ X_1\subset X_0\subset X_{-1} $ and $ X_1 $ is compactly embedded in $ X_0 $. Compactness of $ \mathcal{K} $ in $ L^p(0,T;X_0) $ follows if the set $ \mathcal{K} $ is bounded in $ L^p(0,T;X_1) $ and, as $ \Delta_t\to 0 $, $ \left\|u(\cdot+\Delta_t)-u\right\|_{L^p(0,T-\Delta_t;X_{-1})}\to 0 $, uniformly for $ u\in \mathcal{K} $ [37,Theorem 5].

2.2.   Two-scale convergence

Recall that $ \Omega\subset \mathbb{R} ^3 $ denotes the entire (connected, bounded, open) cardiac domain, assumed to be of the form (1.2). The assumption (1.2) simplifies the presentation. With mild modifications of the upcoming proofs, the results remain valid for general domains with Lipschitz boundary. Let $ Y $ be a reference unit cell in $ \mathbb{R} ^3 $, which we fix to be the unit cube $ Y : = [0,1]^3 $.

Let $ Y_i $ and $ Y_e $ be the (disjoint, connected, open) intra and extracellular spaces within $ Y $, separated by the cell membrane $ \Gamma $:

Denote by $ K^{ \varepsilon } $ the set of $ k\in \mathbb{Z} ^3 $ for which $ \bigcup_{k\in K^ \varepsilon } \varepsilon \left(k+Y \right) = \overline{\Omega} $. We define the intracellular domain $ \Omega_i^{ \varepsilon } $, the extracellular domain $ \Omega_e^{ \varepsilon } $, and the cell membrane $ \Gamma^{\varepsilon} $ as

Both sets $ \Omega_i^{\varepsilon} , \Omega_e^{\varepsilon} $ are connected Lipschitz domains, see Figure 1.1. Note however that it is impossible to have both $ \Omega_i^{\varepsilon} $ and $ \Omega_e^{\varepsilon} $ connected in a two-dimensional picture.

To derive estimates for the microscopic model, we employ the following trace inequality for $ \varepsilon $-periodic hypersurfaces:

for some constant $ C $ independent of $ \varepsilon $, cf. [24,Lemma 3] or [29,Lemma 4.2].

We need a uniform Poincaré inequality for perforated domains [10].

Lemma 2.1. There exists a constant $ C $, independent of $ \varepsilon >0 $, such that

for all $ u \in H^1\left( \Omega_j^{\varepsilon} \right) $, $ 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 $ \left\{u^{ \varepsilon }\right\}_{ \varepsilon >0} \subset L^2((0,T)\times \Omega) $ two-scale converges to $ u $ in $ L^2((0,T)\times \Omega ;L^2_{{\rm{per}}}(Y)) $ if

for all $ \varphi \in C([0,T]\times \overline{\Omega};C_{{\rm{per}}}(Y)) $. We express this symbolically as

By density properties, the convergence (2.10) also holds for test functions $ \varphi $ from $ L^2((0,T)\times \Omega;C_{{\rm{per}}}(Y)) $ [12,p. 176].

The two-scale compactness theorem [1,12] is of fundamental importance.

Theorem 2.2 (two-scale compactness). Let $ \{u^{ \varepsilon }\}_{ \varepsilon >0} $ be a bounded sequence in $ L^2((0,T)\times \Omega ) $, $ \left\|u^{ \varepsilon }\right\|_{L^2((0,T)\times \Omega )} \leq C $ $ \forall \varepsilon >0 $. Then there exist a subsequence $ \varepsilon _n\to 0 $ and a function $ u\in L^2( (0,T)\times \Omega ;L^2(Y)) $ such that $ u^{ \varepsilon _n} $ two-scale converges to $ u $ as $ n \rightarrow \infty $.

Consider a sequence $ \{u_j^{ \varepsilon }\}_{ \varepsilon >0} $ of functions defined on the perforated domain $ (0,T)\times \Omega_j^{\varepsilon} $, $ j = i,e $. We write $ \widetilde{u_j^ \varepsilon } $ for the zero-extension of $ u_j^ \varepsilon $ to $ (0,T)\times \Omega $:

By Theorem $ 2.2 $, $ \left\{\widetilde{u_j^ \varepsilon }\right\}_{ \varepsilon >0} $ has a two-scale convergent subsequence, provided we know that $ \|u_j^{ \varepsilon }\|_{L^2((0,T)\times \Omega_j^{\varepsilon} )}\leq C $ $ \forall \varepsilon >0 $. However, this is not true in general for the gradient of $ \widetilde{u_j^{ \varepsilon }} $, even if $ \left\|u_j^ \varepsilon \right\|_{L^2(0,T;H^1( \Omega_j^{\varepsilon} ))}\leq C $, since the extension by zero creates a discontinuity across $ \Gamma^{\varepsilon} $. Instead the following statement holds true for the gradient:

Lemma 2.3. Fix $ j\in\{i,e\} $ and suppose $ u^{ \varepsilon } = $ $ u_j^ \varepsilon $ satisfies $ \left\|u^{ \varepsilon }\right\|_{L^2(0,T;H^1( \Omega_j^{\varepsilon} ))} \leq C $ $ \forall \varepsilon >0 $. Then there exist a subsequence $ \varepsilon _n\to 0 $ and functions $ u\in L^2(0,T;H^1( \Omega )) $, $ u_1\in L^2((0,T)\times \Omega ;H_{{\rm{per}}}^1(Y_j)) $ such that as $ n\to \infty $,

Here, $ {1}_{Y_j}(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 $ \Gamma^{\varepsilon} $ is given by

where $ \Gamma \subset Y $ is a surface in the unit cell. Since $ | \Gamma^{\varepsilon} | \sim \varepsilon ^{-1} $, it is necessary to introduce a normalizing factor in the definition of two-scale convergence on surfaces.

A sequence $ \{v^{ \varepsilon }\}_{ \varepsilon >0} $ of functions in $ L^2((0,T)\times \Gamma^{\varepsilon} ) $ two-scale converges to $ v $ in $ L^2((0,T)\times \Omega ; L^2(\Gamma)) $, written

if, for all $ \varphi \in C([0,T]\times \overline{ \Omega };C_{{\rm{per}}}( \Gamma )) $,

As with (2.10), this convergence continues to hold for test functions $ \varphi $ that belong to $ L^2((0,T)\times\Omega;C_{{\rm{per}}}( \Gamma )) $.

There is a version of Theorem 2.2 for functions on periodic surfaces [2].

Theorem 2.4 (two-scale compactness on surfaces). Suppose $ \{v^{ \varepsilon }\}_{ \varepsilon >0} $ is a sequence of functions in $ L^2((0,T)\times \Gamma_ \varepsilon ) $ satisfying

for some function $ C $ that is independent of $ \varepsilon >0 $. Then there exist a subsequence $ \varepsilon _n\to 0 $ and a function $ v \in L^2((0,T)\times \Omega ;L^2( \Gamma )) $ such that as $ n\to \infty $,

One can characterize the two-scale limit of traces of bounded sequences in $ L^2(0,T;H^1( \Omega ^{ \varepsilon })) $ as the trace of the two-scale limit [2].

Lemma 2.5. Fix $ j\in\{i,e\} $. Suppose $ u^{ \varepsilon } = u_j^ \varepsilon $ satisfies $ \left\|u^{ \varepsilon }\right\|_{L^2(0,T;H^1( \Omega_j^{\varepsilon} ))} \leq C $ $ \forall \varepsilon >0 $ and, cf. Lemma 2.3, $ \widetilde{u^{ \varepsilon }} \overset{2}{\rightharpoonup} {1}_{Y_j}(y)u(t,x) $ in $ L^2((0,T)\times \Omega ;L^2(Y)) $. Let

be the trace of $ u^{ \varepsilon } $ on $ \Gamma^{\varepsilon} = \partial \Omega_j^{\varepsilon} \setminus \partial \Omega $. Then, up to a subsequence,

Remark 2.6 In view of Lemma 2.5, we have (in the sense of measures)

2.3.   Unfolding operators

An alternative approach to studying convergence on oscillating surfaces $ \Gamma^{\varepsilon} $ is provided by the boundary unfolding operator [10]. For any $ x\in \mathbb{R} ^3 $, we have the decomposition $ x = \lfloor x\rfloor + \{ x\} $, where $ \lfloor x \rfloor\in \mathbb{Z} ^3 $ and $ \{x\}\in [0,1]^3 $ denotes the integer and fractional parts of $ x $, respectively. For later use, note the following simple properties, which hold for any $ x,\bar x\in \mathbb{R} ^3 $, $ n\in \mathbb{Z} ^3 $:

Applying the above decomposition to $ x/ \varepsilon $ gives

where $ \lfloor x/ \varepsilon \rfloor\in \mathbb{Z} ^3 $, $ \{ x / \varepsilon \}\in [0,1]^3 $.

The boundary unfolding operator $ \mathcal{T} _{ \varepsilon }^b $ is defined by

The advantage of the unfolding operator is that we can formulate questions of convergence in a fixed space $ L^2(0,T;L^2( \Omega \times \Gamma )) $. 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\in [1,\infty) $. We refer to [10] for the definition of the boundary unfolding operator and proofs of the properties listed next.

The boundary unfolding operator $ \mathcal{T} _{ \varepsilon }^b $ is bounded, linear, and satisfies

For any $ Y $-periodic function $ \psi\in L^2( \Gamma ) $, set $ \psi_ \varepsilon (x): = \psi(x/ \varepsilon ) $. Then

For $ v\in L^2(0,T;L^2( \Gamma^{\varepsilon} )) $, we have the integration formula

for a.e. $ t\in (0,T) $, thereby converting an integral over the oscillating set $ \Gamma^{\varepsilon} $ to an integral over the fixed set $ \Omega \times \Gamma $. For $ v\in L^2(0,T;L^2( \Gamma^{\varepsilon} )) $,

for a.e. $ t\in (0,T) $. For any $ v\in L^2(0,T;L^2( \Gamma^{\varepsilon} )) $,

for a.e. $ t\in (0,T) $, and also in $ L^2(0,T;L^2( \Omega \times \Gamma )) $. Suppose $ \{v^{ \varepsilon }\}_{ \varepsilon >0} $ is a sequence of functions in $ L^2((0,T)\times \Gamma^{\varepsilon} ) $ satisfying (2.11). Then

We need also the unfolding operators linked to the domains $ \Omega_i^{\varepsilon} , \Omega_e^{\varepsilon} $ [10]:

The unfolding operator $ \mathcal{T} ^j_{ \varepsilon } $ maps functions defined on the oscillating set $ (0,T)\times \Omega_j^{\varepsilon} $ into functions defined on the fixed domain $ (0,T)\times \Omega \times Y_j $. The operator $ \mathcal{T} _{ \varepsilon }^j $ is bounded, linear and satisfies

For any $ Y $-periodic function $ \psi\in L^2(Y_j) $, set $ \psi_ \varepsilon (x): = \psi(x/ \varepsilon ) $. Then

For $ u \in L^2(0,T;L^2( \Omega_j^{\varepsilon} )) $, we have the integration formula

for a.e. $ t\in (0,T) $. The integration formula implies

for a.e. $ t\in (0,T) $. Let $ u \in L^2(0,T;H^1( \Omega_j^{\varepsilon} )) $. Then

and hence $ \mathcal{T} _{ \varepsilon }^j (u)\in L^2( \Omega ;H^1(Y_j)) $, for a.e. $ t\in (0,T) $:

The unfolding operators $ \mathcal{T} _{ \varepsilon }^b $ and $ \mathcal{T} _{ \varepsilon }^j $ are related in the following sense:

for a.e. $ t\in (0,T) $. Combining (2.20), (2.18), (2.19), and the trace inequality (2.1) in $ H^1( \Omega_j^{\varepsilon} ) $, we obtain

for a.e. $ t\in (0,T) $, where the constant $ C $ is independent of $ \varepsilon $ and $ t $. Whenever it is convenient, we will write $ \mathcal{T} _{ \varepsilon }^b (u) $ instead of $ \mathcal{T} _{ \varepsilon }^b(u|_{ \Gamma^{\varepsilon} }) $.

Next, we consider the local average (mean in the cells) operator

and the piecewise linear interpolation operator [10,Definition 2.5]

Given the Lipschitz regularity of $ Y_j $, the interpolation operator $ Q_ \varepsilon ^j $ satisfies the following estimates [10,Propositions 2.7 and 2.8]:

where $ C $ is a constant that is independent of $ \varepsilon $.

3.   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 "$ \varepsilon $-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\in \mathbb{R} $,

for some constants $ c_I > 0 $ and $ c_{I,1}, c_{I,2}, c_{I,3}, c_{I,4}, c_{H,2}\ge 0 $.

The classical FitzHugh-Nagumo model corresponds to

where $ a\in (0,1) $ and $ k,\epsilon >0 $ are constants.

Repeated applications of Cauchy's inequality yields

for some constants $ \gamma>0 $ and $ \beta\ge 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 $ \frac12 (A+A^\top) $ and its skew-symmetric part $ \frac12 (A-A^\top) $. Recall that in a quadratic form $ z\mapsto Az\cdot z $ the skew-symmetric part does not contribute. Therefore, letting $ \lambda_{\text{min}} $ and $ \lambda_{\text{max}} $ denote respectively the minimum and maximum eigenvalues of the symmetric part of $ A $, we have

For $ \mu>0 $, consider the function $ F^\mu: \mathbb{R} ^2\to \mathbb{R} ^2 $ defined by

Denote by $ \lambda^\mu_{\text{min}}(z) $, $ \lambda^\mu_{\text{max}}(z) $ the minimum, maximum eigenvalues of the symmetric part of the matrix $ \nabla F^\mu(z) $. To ensure that weak solutions are unique, we need an additional assumption on $ I, H $ expressed via $ F^\mu $ [8,p. 479]: $ \exists \mu,\lambda>0 $ such that

One can verify that the FitzHugh-Nagumo model (3.1) obeys (3.3) (with $ \mu = \varepsilon k $).

A consequence of (3.3) is that $ \nabla F^\mu(\bar z) z \cdot z \ge \lambda \left|z\right|^2 $ for all $ \bar z, z\in \mathbb{R} ^2 $. Therefore, writing

it follows that

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 $ \varepsilon $-dependence of the conductivities $ \sigma_j^ \varepsilon $ ($ j = i,e $), the applied currents $ s^{ \varepsilon }_j $ ($ j = i,e $), and the initial data $ v_0^{ \varepsilon }, w_0^{ \varepsilon } $ decouples into a "fast" and a "slow" variable:

for some fixed functions $ \sigma_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 $ \alpha >0 $ such that

for $ j = i,e $. Finally, we assume that each $ \sigma_j $ is symmetric: $ \sigma_j^T = \sigma_j $.

The regularity assumption (35) implies that $ \sigma_j $ is an admissible test function for two-scale convergence [1], which means that

for every two-scale convergent sequence $ \{ \varphi ^{ \varepsilon }\}_{ \varepsilon >0} $, $ \varphi ^{ \varepsilon }\overset{2}{\rightharpoonup} \varphi $. This convergence still holds if the second part of (3.5) is replaced by $ \sigma_j \in L^2\left(\Omega; C_{{\rm{per}}}(Y_j)\right) $.

For the stimulation currents we assume the compatibility condition

and the boundedness in $ L^2 $:

which, in view of (3.4), is guaranteed if we take [12,p. 174]

Similarly, we assume that

Throughout this paper we denote by $ C $ a generic constant, not depending on the parameter $ \varepsilon $. The actual value of $ C $ may change from one line to the next.

3.2.   Weak solutions

Testing (1.1) against appropriate functions $ \varphi _i, \varphi _e, \varphi _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 $ \partial \Omega $ 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^{ \varepsilon },u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon }) $ of functions satisfying the following conditions:

ⅰ (algebraic relation).

ⅱ (regularity).

ⅲ (initial conditions).

ⅳ (differential equations).

for all $ \varphi _j\in L^2(0,T;H^1( \Omega_j^{\varepsilon} )) $ with $ \varphi _j\in L^4((0,T)\times \Gamma^{\varepsilon} ) $ ($ j = i,e $), and for all $ \varphi _w\in L^2(0,T;L^2( \Gamma^{\varepsilon} )) $.

Remark 3.3. In (3.13), (3.14) we use $ \left\langle \cdot, \cdot \right \rangle $ to denote the duality pairing between $ H^{-1/2}( \Gamma^{\varepsilon} )+L^{4/3}( \Gamma^{\varepsilon} ) $ and $ H^{1/2}( \Gamma^{\varepsilon} )\cap L^4( \Gamma^{\varepsilon} ) $. For a motivation of the regularity conditions in Definition 3.2, see Remark 3.9 below.

By the embedding (2.6), $ v^{ \varepsilon },w^{ \varepsilon }\in C\left(0,T;L^2( \Gamma^{\varepsilon} ) \right) $ and therefore the pointwise evaluations $ v(0) $, $ w(0) $ in (3.12) are well defined. The time derivative $ \partial_t v^ \varepsilon $ is a distribution belonging to $ L^2(0,T;H^{-1/2}( \Gamma^{\varepsilon} ))+L^{4/3}((0,T)\times \Gamma^{\varepsilon} ) $ with initial values $ v^ \varepsilon _{0} $, so that the integration-by-parts formula (2.4) holds. Consequently, we may replace

in (3.13) and (3.14) by

for all test functions $ \varphi _j \in C^\infty_0([0,T)\times \Omega_j^{\varepsilon} ) $, or for all $ \varphi _j\in L^2(0,T;H^1( \Omega_j^{\varepsilon} )) $ such that $ \,{\partial _t} \varphi _j\in L^2((0,T)\times \Gamma^{\varepsilon} ) $, $ \varphi _j\in L^4((0,T)\times \Gamma^{\varepsilon} ) $, and $ \varphi _j(T) = 0 $. Later, when passing to the limit $ \varepsilon \to 0 $ in (3.13) and (3.14), we make use of the form (3.16).

Remark 3.4. Consider a weak solution $ (u_i^ \varepsilon ,u_e^ \varepsilon ,v^{ \varepsilon },w^{ \varepsilon }) $ according to Definition 3.2. Thanks to (2.1), the trace of $ \varphi _j \in L^2(0,T;H^1( \Omega_j^{\varepsilon} )) $ belongs to $ L^2(0,T;H^{1/2}( \Gamma^{\varepsilon} )) $. By the Sobolev inequality (2.3), the trace of $ \varphi _j $ belongs also to $ L^2(0,T;L^4( \Gamma^{\varepsilon} )) $. In Definition 3.2 we ask additionally that the trace of $ \varphi _j $ belongs to $ L^4((0,T)\times \Gamma^{\varepsilon} ) $ to ensure that the surface terms in (3.13), (3.14) are well-defined

Indeed, $ \int_0^T \left\langle \,{\partial _t}v^{ \varepsilon }, \varphi _j\right\rangle \, dt $ is well-defined for such $ \varphi _j $. Moreover,

For the membrane model (GFHN), the growth condition on $ I $ implies

and therefore $ \left\|I(v^{ \varepsilon },w^{ \varepsilon })\right\|_{L^{4/3}((0,T)\times \Gamma^{\varepsilon} )}^{4/3}<\infty $. Consequently, $ J<\infty $.

We actually have a more precise bound. As $ H^{1/2}( \Gamma^{\varepsilon} )\subset L^4( \Gamma^{\varepsilon} ) $, we have \newline $ H^{-1/2}( \Gamma^{\varepsilon} )\subset L^{4/3}( \Gamma^{\varepsilon} ) $ and

Integrating this over $ t\in (0,T) $ yields $ I(v^{ \varepsilon },w^{ \varepsilon })\in L^{4/3}(0,T;H^{-1/2}( \Gamma^{\varepsilon} )) $.

The integral on the right-hand side of (3.15) can be treated similarly, since $ \left|H(v,w)\right|^2 \le C\left(\left|v\right|^4+\left|w\right|^2\right) $. 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.

Theorem 3.5 (existence of weak solution for microscopic system). Fix any $ \varepsilon >0 $, and suppose (GFHN), (1.2), (3.4), (3.5), (3.6), (3.7), (3.9), and (3.10) hold. Then the microscopic bidomain model (1.1), (1.3) possesses a weak solution $ (u_i^{ \varepsilon },u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon }) $ (in the sense of Definition 3.2). This weak solution satisfies the a priori estimates collected in Lemma 3.8 below.

Remark 3.6. The unknowns $ u_i^ \varepsilon ,u_e^ \varepsilon $ are determined uniquely up to a constant. We can fix this constant by imposing the normalization condition $ \int u_e^ \varepsilon \, dx = 0 $, as is done in Definition 3.2. Weak solutions are unique (and $ L^2 $ stable) provided we impose the structural condition (3.3) on the membrane functions $ I,H $.

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 $ \varepsilon $-periodic surfaces (2.8). A similar result is used in [34].

Lemma 3.7. Let $ u_j\in H^1( \Omega_j^{\varepsilon} ) $, $ j = i,e $, $ \int_{ \Omega_e^{\varepsilon} } u_e \,dx = 0 $, and set $ v: = u_i\big|_{ \Gamma^{\varepsilon} }-u_e\big|_{ \Gamma^{\varepsilon} } $. There is a positive constants $ C $, independent of $ \varepsilon $, such that

Proof. First, since $ \int u_e = 0 $,

To estimate the $ L^2 $ norm of $ u_i $, write $ u_i = \bar{u}_i+\tilde{u}_i $, where $ \bar{u}_i: = \frac{1}{\left| \Omega_i^{\varepsilon} \right|}\int_{ \Omega_i^{\varepsilon} } u_i\, dx $ is constant in $ \Omega_i^{\varepsilon} $ and $ \tilde{u}_i: = u_i-\bar{u}_i $ has zero mean in $ \Omega_i^{\varepsilon} $. Clearly,

In view of the Poincaré inequality (2.9),

Let us bound $ \left\|\bar{u}_i\right\|_{L^2( \Omega_i^{\varepsilon} )}^2 = \frac{\left| \Omega_i^{\varepsilon} \right|}{\left| \Gamma^{\varepsilon} \right|} \left\|\bar{u}_i\right\|_{L^2( \Gamma^{\varepsilon} )}^2 $. Since $ \left| \Gamma^{\varepsilon} \right| = \left|K^{ \varepsilon }\right| \int_{ \varepsilon \Gamma } \,dS = \varepsilon ^{-3} \varepsilon ^2 \int_{ \Gamma }\,dS $ (recall (2.7)), $ \varepsilon \left| \Gamma^{\varepsilon} \right| \geq c >0 $. Because of this and $ \left| \Omega_i^{\varepsilon} \right|\le \left| \Omega \right| $,

Noting that

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 ($ \varepsilon $-independent) a priori estimates.

Lemma 3.8 (basic estimates for microscopic system). Referring Theorem 3.5, the weak solution $ (u_i^{ \varepsilon },u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon }) $ of (1.1), (1.3) satisfies the a priori estimates

where $ C $ is a positive constant independent of $ \varepsilon $.

Proof. We only outline a proof of these (mostly standard) estimates.

Specifying $ ( \varphi _i, \varphi _e, \varphi _w) = (u_i^{ \varepsilon },-u_e^{ \varepsilon },w^{ \varepsilon }) $ as test functions in (3.13), (3.14), and (3.15), adding the resulting equations, applying the chain rule (2.5), and using (3.6), we obtain

Thanks to (3.2),

for some $ \varepsilon $-independent constants $ \gamma>0 $ and $ \beta\ge 0 $.

Using Cauchy's inequality ("with $ \delta $"), the source term can be bounded as

with $ \delta>0 $ small and $ C_\delta $ independent of $ \varepsilon $. Lemma 3.7 and (3.18) ensure that

Insert (3.22), (3.23), (3.24) into (3.21), use (3.6), and choose $ \delta $ (appropriately) small. Integrating the result over the time intervall $ (0,\tau) $, for $ \tau\in (0,T) $, yields

for some positive constant $ C_2 $ independent of $ \varepsilon $. Applying Gröwall's inequality, recalling (3.10) and (3.8), we obtain estimates (a), (c), (d), (e) in (3.20). Estimate (b) follows from (3.24) and (a), (c). Finally, note that (GFHN) implies the bound

Estimate (f) follows by integrating this over $ (0,T)\times \Gamma^{\varepsilon} $ and using (d), (e).

Remark 3.9. Let us motivate the regularity requirements in Definition 3.2 that are not covered by Lemma 3.8. First, due to (3.11), (2.1) and (3.20),

where the constants $ \tilde C_ \varepsilon , C_{ \varepsilon } $ may depend on $ \varepsilon $ (via (2.1) with $ \Gamma $ replaced by $ \Gamma^{\varepsilon} $).

Next, we claim that

for some constant that may depend on $ \varepsilon $. To see this use a version of (3.13) or (3.14) (with time-independent test functions) to write

for a.e. $ t\in (0,T) $ and for any $ \psi\in H^{1/2}( \Gamma^{\varepsilon} ) $ with $ \left\|\psi\right\|_{H^{1/2}( \Gamma^{\varepsilon} )} = 1 $, where

and $ \mathcal{I}_j(\cdot) $ is the right inverse of the trace operator relative to $ \Omega_j^{\varepsilon} $, for $ j = i $ or $ j = e $. Clearly, using the Cauchy-Schwarz inequality and (2.2),

where the constants $ \tilde C_1, \tilde C_2 $ may depend on $ \varepsilon $ (via the inverse trace inequality (2.2) with $ \Gamma $ replaced by $ \Gamma^{\varepsilon} $). As a result, for a.e. $ t\in (0,T) $,

Integrating this over $ t\in (0,T) $ and using (3.20)-(a), (3.8) it follows that

for some constants $ C_1, C_2 $ (that may depend on $ \varepsilon $).

As in Remark 3.4,

and hence, for a.e. $ t\in (0,T) $,

Integrating over $ t\in (0,T) $ and using estimates (3.20)-(d, e), we conclude

for some constant $ C_3 $ (independent of $ \varepsilon $).

For the upcoming convergence analysis, we need a temporal translation estimate for the membrane potential $ v^ \varepsilon $.

Lemma 3.10 (temporal translation estimate in $ L^2 $). Referring back to Theorem 3.5, let $ (u_i^{ \varepsilon },u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon }) $ be the weak solution of (1.1), (1.3). Then $ v^ \varepsilon $ satisfies the following $ L^2 $ temporal translation estimate for some $ \varepsilon $-independent constant $ C>0 $, for any sufficiently small $ \Delta_t>0 $:

Proof. The translated functions $ (u_i^{ \varepsilon },u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon })(t+\Delta_t) $ constitute a weak solution of (1.1) on $ (0,T-\Delta_t) $ with initial data $ (v^ \varepsilon ,w^ \varepsilon )(\Delta_t) $ and stimulation currents $ (s^ \varepsilon _i,s^ \varepsilon _e)(t+\Delta_t) $. We subtract the equations (linked to $ \Omega_i^{\varepsilon} $ and $ \Omega_e^{\varepsilon} $) 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 $ \varphi _i, \varphi _e, \varphi _w $ as

we obtain

Let us write this equation as

By the Cauchy-Schwarz and Minkowski integral inequalities,

Similarly, $ \left|J_2\right|\leq C_2\Delta_t $.

We need the following facts involving a real-valued function $ f\in 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,

Therefore,

An application of Hölder's inequality yields

where we have repeated the argument for $ J_1 $ involving (3.26), with the implication that $ \left\|\int_{t}^{t+\Delta_t}v^{ \varepsilon }(s)\,ds\right\|_{L^4_{t,x}} $ is bounded in terms of $ \Delta_t \left\|v^{ \varepsilon }\right\|_{L^4_{t,x}} $. Moreover, we have used the basic inequality $ ab\le \frac34 a^{4/3} + \frac14 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 (2), temporal estimates are not available for the intra- and extracellular potentials $ u_i^ \varepsilon , u_e^ \varepsilon $.

4.   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:

Definition 4.1 (weak formulation of bidomain system). A weak solution to (1.4), (4.1) is a collection $ (u_i,u_e,v,w) $ of functions satisfying the following conditions:

ⅰ (algebraic relation).

ⅱ (regularity).

ⅲ (initial conditions).

ⅳ (differential equations).

for all test functions $ \varphi _i, \varphi _e\in L^2(0,T; H^1( \Omega ))\cap L^4((0,T)\times \Omega ) $, $ \varphi _w\in L^2(0,T;L^2( \Omega )) $. We denote by $ \left\langle \cdot, \cdot \right\rangle $ the duality pairing between $ L^2(0,T; H^{-1}( \Omega )) + L^{4/3}((0,T)\times \Omega ) $ and $ L^2(0,T; H^1( \Omega ))\cap L^4((0,T)\times \Omega ) $.

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].

Theorem 4.2 (well-posedness of bidomain system). Suppose (1.2) holds, $ I $ and $ H $ satisfy the conditions in (GFHN) and (3.3) (for uniqueness), $ M_i(x) $ and $ M_e(x) $ are bounded, positive definite matrices, $ s_i,s_e\in L^2((0,T)\times \Omega ) $, and $ v_0,w_0 \in L^2( \Omega ) $. Then there exists a unique weak solution to the bidomain system (1.4), (4.1).

We are now in a position to state the main result, which should be compared to Theorem 1.3 in [34].

Theorem 4.3 (convergence to the bidomain system). Suppose conditions (GFHN), (1.2), (3.3), (3.4), (3.5), (3.6), (3.7), (3.9), and (3.10) hold. Let $ \varepsilon $ take values in a sequence tending to zero (e.g. $ \varepsilon ^{-1}\in \mathbb{N} $). Then the sequence $ \left\{u_i^{ \varepsilon },u_e^{ \varepsilon },v^{ \varepsilon },w^{ \varepsilon }\right\}_{ \varepsilon >0} $ of weak solutions to the microscopic system (1.1), (1.3) two-scale converges (in the sense of (4.4) below) to the weak solution $ (u_i,u_e,v,w) $ of the macroscopic bidomain system (1.4), (4.1). Moreover, $ \left\{v^{ \varepsilon }\right\}_{ \varepsilon >0} $ converges strongly in the sense that

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

where $ M_j $ is defined in (1.5), $ j = i,e $.

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^{ \varepsilon } $ (defined on $ \Omega _i^ \varepsilon $) and $ u_e^{ \varepsilon } $ (defined on $ \Omega _e^ \varepsilon $) do not converge strongly, some kind of strong compactness is expected for the $ \varepsilon $-scaled version of the transmembrane potential $ v^{ \varepsilon } = u_i^{ \varepsilon }\big|_{ \Gamma^{\varepsilon} }-u_e^{ \varepsilon }\big|_{ \Gamma^{\varepsilon} } $ (defined on $ \Gamma^{\varepsilon} $), 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 $ \mathcal{T} _{ \varepsilon }^b $, cf. (2.12), to transform the problem of convergence on the oscillating set $ \Gamma^{\varepsilon} $ to the fixed set $ \Omega \times \Gamma $. Roughly speaking, (3.20) is used to conclude that $ \mathcal{T} _{ \varepsilon }^b(v^ \varepsilon ) $ is uniformly bounded in $ L^2_tL^2_xH^{1/2}_y $. In addition, in view of (3.25), $ \mathcal{T} _{ \varepsilon }^b(v^ \varepsilon ) $ possesses an $ \varepsilon $-uniform temporal translation estimate with respect to the $ L^2_{t,x,y} $ norm. As a result, the Simon compactness result (cf. Subsection 2.1) implies that $ \left\{(t,y)\mapsto \mathcal{T} _{ \varepsilon }^b(v^ \varepsilon )(t,x,y)\right\}_{ \varepsilon >0} $ is precompact in $ L^2((0,T)\times \Gamma ) $, for fixed $ x $. Next we demonstrate that $ \mathcal{T} _{ \varepsilon }^b(v^ \varepsilon ) $ is equicontinuous in $ x $ (with values in $ L^2((0,T)\times \Gamma ) $). Applying the Simon-type compactness criterion found in [19] (cf. Theorem 4.6 below), it follows that $ \left\{ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })\right\}_{ \varepsilon >0} $ converges along a subsequence. Owing to the uniqueness of solutions to the bidomain system (4.1), the entire sequence $ \left\{ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })\right\}_{ \varepsilon >0} $ converges (not just a subsequence). We refer to Subsection 4.2 for details. For inspirational works deriving macroscopic models by combining two-scale and unfolding techniques, we refer to [19,18,20,29,32].

4.1.   Extracting two-scale limits

Recall that $ \tilde{\cdot} $ denotes the extension to $ \Omega $ by zero, and that $ {1}_{Y_j} $ is the indicator function of $ Y_j $ ($ j = i,e $). Using the a priori estimates provided by Lemma 3.8, we can apply Lemma 2.3, Theorem 2.4, and Lemma 2.5 to extract subsequences (not relabelled) such that

for some limits $ u_i,u_e \in L^2(0,T,H^1( \Omega )) $, $ u_i^1,u_e^1 \in L^2((0,T)\times \Omega ;H_{{\rm{per}}}^1(Y)) $, and $ w \in L^2((0,T)\times \Omega \times \Gamma ) $. Here we identify $ v = u_i-u_e $ as an element in $ L^2(0,T,H^1( \Omega )) $. It is easily verified that the two-scale limit $ u_e $ satisfies $ \int_{ \Omega } u_e \, dx = 0 $.

Nonlinear functions are not continuous with respect to weak convergence, which prevents us from immediately making the identifications

Using the two-scale convergences in (4.4), Remark 3.3, the choice of test function

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:

and

for all $ \varphi \in C^{\infty}_0((0,T)\times \Omega ) $.

In (4.5) and (4.6), $ M_i $ and $ M_e $ are the homogenized conductivity tensors $ (1.5) $. Let us briefly recall how one arrives at the homogenized conductivities. Setting $ \varphi \equiv 0 $ and considering

as a test function for two-scale convergence, we have by (4.4) that

Note that the oscillating $ \varphi _1 $ term is suppressed in the limit of the weak formulation (3.2) by the $ \varepsilon $-factor, except in the term where a gradient hits the test function. Thus, the two-scale limit $ (u_j,u_j^1) $ satisfies the equation

for all $ \varphi _1 \in C_0^{\infty}((0,T)\times \Omega ;C_{{\rm{per}}}^{\infty}(Y)) $. This equation is satisfied by $ u_j^1 = \chi_j \cdot \nabla_x u_j $, where $ \chi_j $ is the first order corrector (1.6). Hence, for any $ y $-independent function $ \varphi _1(t,x,y)\equiv \psi(t,x) \in C_0^{\infty}((0,T)\times \Omega ) $,

so (1.5) is indeed the homogenized conductivity tensor.

Up to the convergences of the nonlinear terms, it is now clear that (4.5), (4.6), and (4.7) constitute the weak formulation of the bidomain system (1.4), (4.1) (in the sense of Definition 4.1).

4.2.   The nonlinear terms and strong convergence

To finalize the proof of Theorem 4.3, it remains to identify the limits

for all $ \varphi \in C^{\infty}_0((0,T)\times \Omega ) $. We note that the unfolding operator (2.14) allows us to transform the oscillating surface integral

into

coming from the integration formula (2.14) for $ \mathcal{T} _{ \varepsilon }^b $. Additionally, we have here used (2.13) and the fact $ \mathcal{T} _{ \varepsilon }^b \left(I(v^{ \varepsilon },w^{ \varepsilon })\right) = I\left( \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon }), \mathcal{T} _{ \varepsilon }^b(w^{ \varepsilon })\right) $. The smoothness of $ \varphi $ implies that $ \mathcal{T} _{ \varepsilon }^b( \varphi )\to \varphi $ in $ L^2((0,T)\times \Omega \times \Gamma ) $ as $ \varepsilon \to 0 $, cf. (2.16), so to identify the limits (4.8) it suffices to show

where $ v,w $ are identified in (4.4). Besides, since $ w^{ \varepsilon } $ appears linearly in $ I $ and $ H $, cf. (GFHN), the weak convergence of $ \left\{ \mathcal{T} _{ \varepsilon }^b(w^{ \varepsilon })\right\}_{ \varepsilon >0} $ in $ L^2((0,T)\times \Omega \times \Gamma ) $ is enough to pass to the limit in (4.8), if we establish strong convergence of $ \left\{ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })\right\}_{ \varepsilon >0} $.

As a step towards verifying the required strong convergence, we need to show that $ \left\{(t,y)\mapsto \mathcal{T} _{ \varepsilon }^b (v^ \varepsilon )\right\}_{ \varepsilon >0} $ is strongly precompact in $ L^2_{t,y} $, for fixed $ x\in \Omega $. As a result of Lemma 3.8, $ \mathcal{T} _{ \varepsilon }^b (v) $ is bounded in $ L^2_tL^2_xH^{1/2}_y $, uniformly in $ \varepsilon $. However, according to Lemma 3.8, the time derivative $ \,{\partial _t}v^ \varepsilon $ is merely of order $ 1/ \varepsilon $ in the $ L^2_tH^{-1/2}_x $ norm. Therefore, we cannot expect $ \,{\partial _t} \mathcal{T} _{ \varepsilon }^b(v^ \varepsilon ) $ (assuming that this object is meaningful) to be bounded in $ L^2_tL^2_xH^{-1/2}_y $, uniformly in $ \varepsilon $. As a consequence, strong $ L^2_{t,y} $ compactness of $ \left\{ \mathcal{T} _{ \varepsilon }^b (v^ \varepsilon )\right\}_{ \varepsilon >0} $ is not deducible from the classical Aubin-Lions theorem. Instead of attempting to control (in a negative space) the whole derivative $ \,{\partial _t} \mathcal{T} _{ \varepsilon }^b(v^ \varepsilon ) $, 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 $ \left\{ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })\right\}_{ \varepsilon >0} $ is strongly precompact in the fixed space $ L^2((0,T)\times \Omega \times \Gamma ) $. We begin with

Lemma 4.5. There exists a constant $ C $, independent of $ \varepsilon $, such that

and

for sufficiently small temporal shifts $ \Delta_t>0 $.

Proof. By (3.11) and the linearity of $ \mathcal{T} _ \varepsilon ^b(\cdot) $, we have $ \mathcal{T} _{ \varepsilon }^b (v^{ \varepsilon }) = \mathcal{T} _{ \varepsilon }^b(u_i^{ \varepsilon }|_{ \Gamma^{\varepsilon} }) - \mathcal{T} _{ \varepsilon }^b(u_e^{ \varepsilon }|_{ \Gamma^{\varepsilon} }) $. Thus, (4.9) follows by integrating (2.21) over $ (0,T) $ and using estimates (a), (b) in (3.20). Regarding (4.10), we use the linearity of $ \mathcal{T} _{ \varepsilon }^b $, (2.15), and (3.25) to obtain

Let us think of $ \mathcal{T} _{ \varepsilon }^b(v^ \varepsilon ) $ as a function of $ x\in \Omega $, with values in $ L^2(0,T;L^2( \Gamma )) $. Fixing $ x $, in view of Lemma 4.5 and Simon's compactness criterion, the sequence $ \left\{(t,y)\mapsto \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })(t,x,y)\right\}_{ \varepsilon >0} $ is precompact in $ L^2((0,T)\times \Gamma ) $. The $ x $-variable is more difficult. As a matter of fact, since $ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon }) $ is piecewise constant as a function of $ x $ and thus does not belong to any Sobolev space, strong compactness in $ x $ is not immediately clear. We address this issue by deriving a translation estimate in the $ x $-variable. To be more precise, we make use of a convenient Simon-type compactness criterion ($ x $ playing the role of time) established in [19,Corollary 2.5] (see also [3,Section 5]), which is recalled next.

For $ Q\subset \mathbb{R} ^n $ and $ \xi\in \mathbb{R} ^n $ ($ n\ge 1 $), we set $ Q_\xi: = Q \cap (Q-\xi): = \left\{x\in Q: x+\xi\in Q\right\} $, $ \Sigma: = \left\{-1,1\right\}^n $ ($ \left|\Sigma\right| = 2^n $), and $ \xi_\sigma: = (\xi_1 \sigma_1,\ldots,\xi_n\sigma_n)\in \mathbb{R} ^n $ for $ \sigma\in \Sigma $. If $ Q = (a,b) $ ($ : = \Pi_{\ell = 1}^n (a_\ell,b_\ell) $ for $ a,b\in \mathbb{R} ^n $, $ a<b $) is an open rectangle, then

Let $ B $ be Banach a space. For $ f:Q\to B $ and $ \Delta\in \mathbb{R} ^n $, we define the translation operator $ \tau_\Delta:(Q-\Delta)\to B $ by

The following theorem, due to Gahn and Neuss-Radu [19], is a multi-dimensional generalization of Simon's main result [37,Theorem 1].

Theorem 4.6 ([19]). Let $ \mathcal{F} \subset L^p(Q;B) $ for some Banach space $ B $, open rectangle $ Q = (a,b)\subset \mathbb{R} ^n $, and $ p\in [1,\infty) $. Then $ \mathcal{F} $ is precompact in $ L^p(Q;B) $ if and only if

i. $ \left\{\int_A f \,dx \, \big | \, f \in \mathcal{F}\right\} $ is precompact in $ B $, for every open rectangle $ A\subset Q $;

ii. for each $ z\in \mathbb{R} ^n $ with $ 0<z<b-a $,

Recall (4.11), this time specifying $ Q = (a,b)\subset \mathbb{R} ^n $, and $ \xi = (b-a)/2\in \mathbb{R} ^n $. Condition (4.12) in Theorem 4.6 is equivalent to [19]

The difference between (4.12) and (4.13) is the fixed domain that is utilized in the latter (it does not depend on the shift $ z $). We make use of (4.13) in the proof of Lemma 4.8 below.

We now verify that the sequence $ \left\{ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })\right\}_{ \varepsilon >0} $ of unfolded membrane potentials satisfies the assumptions of Theorem 4.6, with $ B = L^2((0,T)\times \Gamma ) $, $ Q = \Omega $, $ p = 2 $.

Lemma 4.7 (verification of $ i $). Given an arbitrary open rectangle $ A\subset \Omega $, define the function $ v^{ \varepsilon }_A(t,y) $ by

Then the sequence $ \left\{v^ \varepsilon _A\right\}_{ \varepsilon >0} $ is precompact in $ L^2((0,T)\times \Gamma ) $.

Proof. In view of Jensen's inequality, it follows that

Let $ \Delta_t>0 $ be a small temporal shift. Again using Jensen's inequality,

Summarizing, there exists an $ \varepsilon $-independent constant $ C $ such that

The lemma follows from these estimates and Simon's compactness criterion.

In the next lemma we verify (4.13) with $ \xi = (1/2,1/2,1/2)\in \mathbb{R} ^3 $, which is equivalent to condition $ ii $ in Theorem 4.6.

Lemma 4.8 (verification of $ ii $). Given any $ \delta>0 $, there exists $ h>0 $ such that for any $ \Delta_x\in \mathbb{R} ^3 $ with $ \left|\Delta_x\right|<h $ and for all $ \varepsilon \in (0,1] $,

Proof. We wish to estimate the quantity

Recall that $ \varepsilon $ takes values in a sequence $ \subset (0,1] $ tending to zero. Fix any $ \varepsilon _0>0 $. Since the translation operation is continuous in $ L^2 $, there exists $ h_0 = h_0( \varepsilon _0)>0 $ such that $ J(\Delta_x; \varepsilon )<\delta $ for any $ \left|\Delta_x\right|<h_0 $, for all $ \varepsilon \in [ \varepsilon _0,1] $. The rest of the proof is devoted to arguing that this holds also for $ \varepsilon \in (0, \varepsilon _0) $, thereby proving (4.14).

Choose $ K^ \varepsilon _{\xi_\sigma}\subset \mathbb{Z} ^3 $ such that $ \Omega _{\xi_\sigma} = \text{interior}\left(\bigcup_{k\in K^ \varepsilon _{\xi_\sigma}} \overline{ \varepsilon Y^k}\right), \qquad Y^k : = k+Y $. Then $ J(\Delta_x, \varepsilon ) $ becomes

If $ x\in \varepsilon Y^k $, then $ \left\lfloor \frac{x}{ \varepsilon }\right\rfloor = k $, but we have no useful information about $ \left\lfloor \frac{x+(\Delta_x)_\sigma}{ \varepsilon }\right\rfloor $. To address this issue, we make use of a favorable decomposition of the cells $ \varepsilon Y^k $ proposed in [32] (and also utilized in e.g. [19,18]).

We decompose each cell $ \varepsilon Y^k $ as

for $ k\in K^ \varepsilon _{\xi_\sigma} $ and $ \sigma\in\Sigma $. Regarding the translation, for $ x\in \varepsilon Y^{k,m}_\sigma $, we write $ (\Delta_x)_\sigma = \varepsilon \left\lfloor \frac{(\Delta_x)_\sigma}{ \varepsilon }\right\rfloor + \varepsilon \left\{\frac{(\Delta_x)_\sigma}{ \varepsilon }\right\} $, and note that

As a result of this,

where we have also used $ \int_{ \varepsilon Y^{k,m}_\sigma}\, dx \le \int_{ \varepsilon Y^k} \,dx = \varepsilon ^3 $ to arrive at the final line. Since $ \sum_{k\in K^ \varepsilon _{\xi_\sigma}}\int_{ \varepsilon (k+ \Gamma )} = \int_{( \Gamma^{\varepsilon} )_{\xi_\sigma}} $, we conclude that

where the shift $ z = z(\Delta_x, \varepsilon ,m) $ is $ \varepsilon \left(m_\sigma +\left\lfloor \frac{(\Delta_x)_\sigma}{ \varepsilon }\right\rfloor\right) $, i.e., $ z $ is an integer-multiple of $ \varepsilon $. Note that $ x+z\in ( \Gamma^{\varepsilon} )_{\xi_\sigma} $ whenever $ \Delta_x $ and $ \varepsilon $ are sufficiently small. Recalling the definition (3.11) of $ v^ \varepsilon $, the trace inequality (2.8) implies

where the last term is bounded by a constant times $ \varepsilon ^2 $ because of (3.20)-(a).

It remains to estimate the term on the second line, which will be done utilizing the well-known characterization of Sobolev spaces by means of translation (difference) operators. Recalling the standard proof of this characterization, a problem that arises (due to the geometry of $ \Omega_j^{\varepsilon} $) is that parts of the line segment between $ x $ and $ z $ may leave $ \Omega_j^{\varepsilon} $. To avoid this problem we make use of the interpolation operators (2.22) to obtain functions $ Q_ \varepsilon ^j(u_j^ \varepsilon ) $ defined on the whole of $ \Omega $.

Using the triangle inequality and recalling the estimates (2.23), we obtain

Hence,

where the constants $ C_5,C_6 $ are independent of $ \Delta_x, \varepsilon $. We select the $ \varepsilon _0 $ introduced earlier sufficiently small, such that the first term on the right-hand side is $ <\delta^2/2 $ for all $ \varepsilon < \varepsilon _0 $. We pick $ h_1>0 $ such that the second term is $ <\delta^2/2 $ for all $ \left|\Delta_x\right|<h_1 $ (for any $ \varepsilon \in (0,1] $). Specifying $ h: = \min(h_0,h_1) $, the claim (4.14) now follows.

4.3.   Concluding the proof of Theorem 4.3

Summarizing, we know that

because $ w^{ \varepsilon }\overset{2-{\rm{S}}}{\rightharpoonup} w $, cf. (4.4) and (2.17). Besides, $ w^{ \varepsilon } $ appears linearly in $ I $ and $ H $, cf. (GFHN). We have shown that $ \left\{ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })\right\}_{ \varepsilon >0} $ is strongly precompact in $ L^2((0,T)\times \Omega \times \Gamma ) $. It then follows that $ \mathcal{T} _{ \varepsilon }^b(v^{ \varepsilon })\to v $ in $ L^2((0,T)\times \Omega \times \Gamma ) $ and a.e. in $ (0,T)\times \Omega \times \Gamma $, along a subsequence as $ \varepsilon \to 0 $ (not relabelled), where $ v $ is the two-scale limit of $ \left\{v^{ \varepsilon }\right\}_{ \varepsilon >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 $ \varepsilon $. 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.