Dominique Chapelle
email: dominique.chapelle@inria.fr

Annabelle Collin
email: annabelle.collin@inria.fr

Jean-Fr Éd

Éric Gerbeau

A SURFACE-BASED ELECTROPHYSIOLOGY MODEL RELYING ON ASYMPTOTIC ANALYSIS AND MOTIVATED BY CARDIAC ATRIA MODELING

Keywords: Computational electrophysiology, Asymptotic analysis, Thin domains, Cardiac modeling AMS Subject Classification: 22E46, 53C35, 57S20

published or not. The documents may come

A surface-based electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling

Introduction

Cardiac electrophysiology purports to describe and model chemical and electrical phenomena taking place in the cardiac tissue, and which are responsible for activating the mechanical contraction in the myocytes, namely, the cardiac muscle cells, see e.g. Refs. 16, 15, 19, 3 and references therein. Given the frequent occurrence of pathologies -such as atrial fibrillation and ventricular tachycardia -directly affecting electrophysiology, hence impairing cardiac function, the detailed understanding of the associated underlying mechanisms is particularly important. In this context, there is a tremendous potential for modeling, in particular with models adequately "personalized" to given patients, in order to guide medical decision in selecting the best-adapted therapeutic strategy. [START_REF] Smith | euHeart: personalized and integrated cardiac care using patient-specific cardiovascular modelling[END_REF][START_REF] Sermesant | Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation[END_REF] For instance, personalized models may be used to optimize the procedure called radio-frequency ablation that aims at treating cardiac arrythmias -e.g., atrial fibrillation -by blocking some abnormal conduction pathways.

Computational electrophysiology modeling -when aiming at representing the behavior of the global heart or whole subparts thereof, which of course is crucial when considering conduction phenomena and their incidence on cardiac function -is known to give rise to highly intensive simulations. [START_REF] Colli Franzone | Simulating patterns of excitation, repolarization and action potential duration with cardiac Bidomain and Monodomain models[END_REF][START_REF] Göktepe | Computational modeling of cardiac electrophysiology: A novel finite element approach[END_REF] This is due, in particular, to the substantial refinements required in the meshes in order to accurately capture the conduction waves. As personalized modeling generally involves the solution of inverse problems in order to estimate anatomical and biophysical parameters characterizing a given patient-specific model -in order to enable predictive modeling -based on clinical data, this may lead to prohibitive computational times, indeed. The same difficulty holds as regards optimization loops needed to adapt the therapeutic strategy.

Nevertheless, when considering the thin-ness of cardiac structures such as the atria, there is clearly a case for modeling simplifications based on considering surface-based descriptions -hence, in essence two-dimensional (2D) -instead of 3D models, very much in the spirit of structural models -such as shells -in mechanics. Such mechanical models have been extensively analyzed and provide extremely computationally-effective formulations, with very limited loss of accuracy when the structure is thin. As regards electrophysiology modeling, this path has been scarcely explored. A frequently-argued explanation for that lies in the complexity of the 3D anatomical details that appear to be needed to accurately capture the relevant phenomena. In particular, the fiber architecture of the muscle gives rise to strong anisotropy in the diffusion phenomena, and the fiber directions are known to be very rapidly varying across the thickness of the cardiac walls, including in the atria. [START_REF] Ho | The Importance of Atrial Structure and Fibers[END_REF] This makes naive surface reductions of the models based on isotropy or "average anisotropy" assumptions very inaccurate, and in fact altogether unable to represent some complex conduction patterns. [START_REF] Krueger | Modeling atrial fiber orientation in patient-specific geometries: a semi-automatic rule-based approach[END_REF] The objective in this paper is to propose a surface-based electrophysiology model relying on a detailed asymptotic analysis. This type of derivation -also available for shell structures in mechanics, see e.g. Refs. 6, 4, 5 -provides a "mathematically justified" surface model designed for thin structures. In our case, in order to incorporate the above-mentioned key anatomical details, we include the rapid variations of the fiber direction within our asymptotic analysis assumptions. We point out that some asymptotic derivations have also been performed in earlier works, [START_REF] Poignard | About the transmembrane voltage potential of a biological cell in timeharmonic regime[END_REF] albeit without taking into account the anisotropy and variations thereof.

asymptoticAnaly-sis˙m3as

A surface-based electrophysiology model relying on asymptotic analysis 3 An outline of the paper is as follows. In Sections 2 and 3 we recall some relevant elements of cardiac electrophysiology modeling, and differential geometry, respectively. Next, in Section 4 we formulate an asymptotic framework for an anisotropic diffusion model with anisotropy directions rapidly varying across the thickness, and we establish strong convergence results to a limit model. This allows us to propose a surface-based electrophysiology model in Section 5, with detailed numerical assessments.

Elements of cardiac electrophysiology

The electric wave propagating in the cardiac tissue can be represented by a nonlinear reaction-diffusion partial differential equation (PDE), coupled with an ordinary differential equation (ODE) representing cellular activity. Considering in particular the bidomain model -see for example Refs. 16, 15, 19, 2 -equations can be written in terms of the extracellular potential u e , the transmembrane potential V m = u i -u e , with u i the intracellular potential, and the ionic variable w, as

             A m C m ∂V m ∂t + I ion (V m , w) -div σ i • ∇V m = div σ i • ∇u e + A m I app , in B × (0, T ), div σ i + σ e • ∇u e = -div σ i • ∇V m , in B × (0, T ), ∂ t w + g(V m , w) = 0, in B × (0, T ), (2.1)
where B denotes the 3D domain of interest, and with appropriate boundary conditions

σ i • ∇u e • n = -σ i • ∇V m • n, in ∂B × (0, T ), σ e • ∇u e • n = 0, in ∂B × (0, T ), (2.2)
where A m is a positive constant denoting the ratio of membrane area per unit volume, C m the membrane capacitance per unit surface, I ion (V m , w) a reaction term representing the ionic current across the membrane, and I app a given applied stimulus current. The current I ion can be described by a phenomenological model. In this study, the model proposed by Mitchell and Schaeffer in Ref. 13 is considered. The functions g and I ion are then given by

I ion (V m , w) = - w(V m -V min ) 2 (V max -V m ) τ in (V max -V min ) + V m -V min τ out (V max -V min ) , g(V m , w) =      w τ open - 1 τ open (V max -V min ) 2 if V m < V gate , w τ close if V m ≥ V gate ,
where τ in , τ out , τ open , τ close are given parameters and V min , V max , V min < V gate < V max auxiliary constants.

Cardiac muscle has a fiber architecture. The electrical conductivity is higher along than across the fiber direction. The intracellular and extracellular media are therefore anisotropic. This anisotropy is included in our model defining the conductivity tensors σ i and σ e by σ i,e = σ t i,e g + (σ l i,eσ t i,e ) τ ⊗ τ , where g denotes the 3D metric tensor -the components of which are given by the identity matrix in an orthonormal coordinate system, see e.g. Ref. 4 -the vector τ is of unit length and parallel to the local fiber direction, and σ l i and σ t i are the conductivity coefficients in the intra-cellular medium measured along and across -respectively -the fiber direction, and likewise for σ l e and σ t e in the extra-cellular medium.

The bidomain model can be rewritten in weak form as follows. For all t > 0, find

V m (•, t) ∈ H 1 , u e (•, t) ∈ H 1 and w(•, t) ∈ L ∞ with B u e = 0, such that                  A m B C m ∂V m ∂t + I ion (V m , w) φ + B σ i • ∇V m + ∇u e • ∇φ = A m B I app (V m , w)φ, B ( σ i + σ e ) • ∇u e • ∇ψ + B σ i • ∇V m • ∇ψ = 0, ∂ t w + g(V m , w) = 0, in B,
for all φ, ψ ∈ H 1 such that B ψ = 0. With our motivation of deriving a 2D model valid for thin cardiac structuresthe atrial walls, in particular -and defined over the midsurface of the thin region, we observe that the main difficulty in representing the anisotropy resulting from the preferred conduction direction along the muscle fibers -which may rapidly vary across the thickness -arises from the diffusion term in the equations. Hence, in our mathematical analysis we will focus on this term, with the primary concern of taking into account both major features, namely, anisotropic behavior and heterogeneous distribution across the thickness.

Geometric preliminaries

The midsurface S of the body is described by a mapping φ defined over ω, with ω a domain of R 2 , and with values in the three-dimensional Euclidean space E, see Figure 1. We assume that φ is as regular as needed. For the sake of simplicity in our asymptotic analysis we suppose that the thickness of the body, denoted by d, is constant over S. We introduce L = diam(S) and ε = d L a dimensionless thickness parameter. Note that this parameter is defined to be geometrically intrinsic, namely, it is invariant with respect to changes in the geometric mapping φ (hence, in the coordinate system and its domain ω) that preserve the surface.

Denoting by (ξ 1 , ξ 2 ) the coordinates used in R 2 (hence in ω), we define

a α (ξ 1 , ξ 2 ) = ∂ α φ(ξ 1 , ξ 2 ), α = 1, 2.

asymptoticAnaly-sis˙m3as

A surface-based electrophysiology model relying on asymptotic analysis 5

We suppose that φ is such that the vectors ( a α ) α=1,2 form a basis -called the covariant basis -for the tangential plane to the midsurface S at any point with coordinates in ω. We also introduce the contravariant basis ( a 1 , a 2 ) such that

a α (ξ 1 , ξ 2 ) • a β (ξ 1 , ξ 2 ) = δ β α , α, β = 1, 2, ∀(ξ 1 , ξ 2 ) ∈ ω,
and we denote by

a αβ = a α • a β ,
the covariant components of the surface-attached metric tensor, see e.g. Ref. 4.

Remark 3.1. Introducing the 2 × 2 matrix A = (a αβ ) α,β=1,2 , we have

A -1 = (a αβ ) α,β=1,2 = ( a α • a β ) α,β=1,2 .
We denote by a the quantity a = a 1 × a 2 2 , and the unit vector normal to the tangential plane at any point with coordinates in ω is given by

a 3 = a 1 × a 2 √ a . Defining now Ω = ω × - d 2 , + d 2
, the 3D geometry is described by the mapping

Φ : Ω -→ E (ξ 1 , ξ 2 , ξ 3 ) -→ φ(ξ 1 , ξ 2 ) + ξ 3 a 3 (ξ 1 , ξ 2 )
and the 3D geometric domain of interest is given by B = Φ Ω , see Figure 1. We also introduce the 3D covariant basis

g i (ξ 1 , ξ 2 , ξ 3 ) = ∂ i Φ(ξ 1 , ξ 2 , ξ 3 ), i = 1, 2, 3, ∀(ξ 1 , ξ 2 , ξ 3 ) ∈ Ω,
and the corresponding contravariant basis ( g 1 , g 2 , g 3 ) The components of the 3D metric tensor g in covariant form are then given by

g i (ξ 1 , ξ 2 , ξ 3 ) • g j (ξ 1 , ξ 2 , ξ 3 ) = δ j i , i, j = 1, 2, 3, ∀(ξ 1 , ξ 2 , ξ 3 ) ∈ Ω.
g ij = g i • g j . Remark 3.2. Introducing the 3 × 3 matrix G = (g ij ) i,j=1,2,3 , we have G -1 = (g ij ) i,j=1,2,3 = ( g i • g j ) i,j=1,2,3 .
For a differentiable scalar field defined over the surface the 2D gradient is given by

∇u = ∂ 1 u a 1 + ∂ 2 u a 2 ,
whereas for a volume-defined quantity we have the 3D gradient

∇u = ∂ 1 u g 1 + ∂ 2 u g 2 + ∂ 3 u g 3 .
We need a number of preliminary results that relate 3D quantities to surface quantities. For the proofs of the three following propositions, see e.g. Ref.

b αβ = a 3 • ∂ β a α , c αβ = λ,µ=1,2 a λµ b λα b µβ ,
we have

g αβ = a αβ -2ξ 3 b αβ + (ξ 3 ) 2 c αβ , (3.1)
and, moreover,

g α3 = 0, g 33 = 1. (3.2)
Proposition 3.2. For the volume measure given by dV = √ g dξ 1 dξ 2 dξ 3 , with g = (det( g 1 , g 2 , g 3 )) 2 = det(G), we have

g = a 1 -2Hξ 3 + K(ξ 3 ) 2 2 ,
denoting by H and K the mean and Gaussian curvatures of the surface S, respectively. Assuming that 1 -2Hξ 3 + K(ξ 3 ) 2 is strictly positive over B, we then infer

√ g = √ a 1 -2Hξ 3 + K(ξ 3 ) 2 and ∃γ > 0, √ g > γ √ a.
Proposition 3.3. There exist two strictly positive constants c 1 , c 2 such that,

∀(ξ 1 , ξ 2 , ξ 3 ) ∈ Ω, ∀(x 1 , x 2 ) ∈ R 2 , c 1 2 α,β=1 a αβ (ξ 1 , ξ 2 )x α x β ≤ 2 α,β=1 g αβ (ξ 1 , ξ 2 , ξ 3 )x α x β ≤ c 2 2 α,β=1 a αβ (ξ 1 , ξ 2 )x α x β .
We can then show the following identity.

Lemma 3.1. ∀α, β = 1, 2, there exist d αβ and e αβ , two functions of (ξ 1 , ξ 2 ) in L ∞ (ω), such that

g αβ = 1 1 -2Hξ 3 + K(ξ 3 ) 2 2 a αβ -2ξ 3 d αβ + (ξ 3 ) 2 e αβ .
April We introduce the space H 1 (B) defined by

H 1 (B) = u : Ω → R measurable, Ω u 2 dV, Ω ∇u • ∇u dV < +∞ ,
with the natural norm

u H 1 (B) = Ω u 2 dV + Ω ∇u • ∇u dV 1 2
.

We need to introduce boundary conditions associated with the diffusion term that we want to study. We suppose ∂ω = γ 1 ∪ γ 2 with γ 1 of non-zero measure, and we define

Γ i = γ i × -d 2 , + d 2 , for i = 1, 2. We denote Γ 3 = ∂Ω \ Γ 1 ∪ Γ 2 ,
namely, the top and bottom surfaces. We consider the following problem

         -div σ • ∇u = f, in Ω, u = 0, in Γ 1 , σ • ∇u • n = 0, in Γ 2 , σ • ∇u • n = 0, in Γ 3 , (4.1)
where f ∈ L 2 (B).

Defining

V 3D = H 1 (B)∩(BC)
, where (BC) corresponds to the Dirichlet boundary condition on Γ 1 , the problem (4.1) can be rewritten in weak form as seeking u ∈ V 3D such that

A 3D (u, v) = F 3D (v), ∀v ∈ V 3D , (4.2)
where

A 3D (u, v) = Ω σ • ∇u • ∇v dV, F 3D (v) = Ω f v dV, ∀u, v ∈ V 3D .
Remark 4.1. The variational problem is here written in an intrinsic form, namely, independently of any specific coordinate system, but of course it is very straightforward to obtain the corresponding expressions using the tensor components, which we will consider soon with the coordinates introduced in the previous section.

Assumption 4.1. We suppose that the diffusion tensor σ(ξ) is:

• L ∞ on Ω;

• symmetric positive definite for almost all ξ in Ω;

• such that for almost all ξ ∈ Ω, for any ρ ∈ E,

c 1 ρ 2 ≤ σ(ξ) • ρ • ρ,
for some strictly positive constant c 1 .

Remark 4.2. Due to the L ∞ character of σ(ξ) we also have

σ(ξ) • ρ • ρ ≤ c 2 ρ 2 ,
for some strictly positive c 2 , for any ρ ∈ E and for almost all ξ ∈ Ω.

A simple application of the Lax-Milgram theorem with the Poincaré inequality then provides the following result.

Theorem 4.1. If Assumption 4.1 holds, there exists a unique u ∈ V 3D solution of (4.2).

Asymptotic problem formulation

Denoting dS = a(ξ 1 , ξ 2 ) dξ 1 dξ 2 the surface measure, we introduce the spaces L 2 (S) and H 1 (S) by

L 2 (S) = u : ω → R measurable, ω u 2 dS < +∞ , H 1 (S) = u : ω → R measurable, ω u 2 dS, ω ∇u • ∇u dS < +∞ , with the natural norms u L 2 (S) = ω u 2 dS 1 2 , u H 1 (S) = ω u 2 dS + ω ∇u • ∇u dS 1 2
.

We assume that the source term f is smooth enough to provide

f (ξ 1 , ξ 2 , ξ 3 ) = f 0 (ξ 1 , ξ 2 ) + ξ 3 f (ξ 1 , ξ 2 , ξ 3 ),
where

f 0 ∈ L 2 (S) and f ∈ L ∞ (B).
Regarding the modeling of the anisotropic conductivity, we make the following assumption on the diffusion tensor σ.

Assumption 4.2. We suppose that σ = σ t g + (σ lσ t ) τ ⊗ τ , where σ t and σ l are two strictly positive constants such that σ t ≤ σ l . The vector τ is given by

τ (ξ 1 , ξ 2 , ξ 3 ) = τ 0 (ξ 1 , ξ 2 ) cos 2θ(ξ 1 , ξ 2 )ξ 3 d + τ ⊥ 0 (ξ 1 , ξ 2 ) sin 2θ(ξ 1 , ξ 2 )ξ 3 d , (4.3)
where τ 0 and τ ⊥

• τ 0 belongs to the tangential plane to the midsurface S, at each point in ω,

• τ ⊥ 0 = a 3 × τ 0 i.e τ ⊥ 0 is orthogonal to τ 0 and belongs to the tangential plane.

Remark 4.3. The assumptions on σ translate the anisotropic behavior and heterogeneous distribution across the thickness of the body, see Figure 2. Note that

• σ l = σ t gives the homogeneous case, namely, without any privileged direction;

• σ t < σ l implies (a) if θ = 0, then τ is independent of ξ 3 , i.e
. the fiber orientation is constant across the thickness, defined by the vector τ 0 ; (b) if θ = 0, the fibers rotate across the thickness by a total angle 2θ.

Remark 4.4. We commit a slight abuse of notation in (4.3), since we have a 3D vector in the left-hand side and surface-attached 2D vectors in the right-hand side. We could have instead used a 3D notation for τ 0 and τ ⊥ 0 , but we made this choice to emphasize that they lie in the tangential plane. An asymptotic analysis then consists in studying the behavior of the above diffusion model when varying the thickness parameter ε while keeping the midsurface S and the parameter field θ fixed. The complete asymptotic analysis of this 3D anisotropic diffusion problem is out of the scope of our preliminary analysis, which primarily aims at guiding the formulation of a relevant surface-based model. Hence, following similar strategies used in shell mechanical modeling, we will instead consider a Galerkin reduction of the problem (4.2) posed in a subspace of V 3D given by polynomial variations of the quantity of interest in the thickness variable -in our τ 0 d -θ +θ Fig. 2. Fibers rotation across the thickness case quadratic polynomials -which can be seen as truncated Taylor expansions. We thus define V = (H 1 (S) ∩ (BC )) [START_REF] Carusi | Bridging experiments, models and simulations: an integrative approach to validation in computational cardiac electrophysiology[END_REF] , where (BC ) corresponds to homogeneous Dirichlet boundary conditions prescribed on γ 1 , and we seek (

u ε 0 , u ε 1 , u ε 2 ) ∈ V such that A 3D u ε 0 + ξ 3 u ε 1 + (ξ 3 ) 2 u ε 2 , v 0 + ξ 3 v 1 + (ξ 3 ) 2 v 2 = F 3D v 0 + ξ 3 v 1 + (ξ 3 ) 2 v 2 , ∀(v 0 , v 1 , v 2 ) ∈ V, (4.4)
where we point out that the definitions of A 3D and F 3D include a dependence on the varying thickness parameter ε -via the dependence of the reference domain Ω itself -even though this is not reflected in the notation for compactness purposes.

Remark 4.6. A 3D asymptotic analysis -namely, of a sequence of solutions of (4.2) associated with decreasing thickness parameter for the geometric domainwould require different tools, starting with a scaling of the domain in the transverse direction in order to work on a fixed domain, see e.g. Ref. 6. Here, we adopt a somewhat simpler strategy with a Galerkin projection of the 3D model onto a subspace of functions that feature quadratic variations in the transverse direction. Note that this is also a classical asymptotic approach in structural mechanics and associated mathematical analysis, where it is known to provide important insight into the limit asymptotic behavior. [START_REF] Delfour | Intrinsic P (2, 1) thin shell model and Naghdi's models without a priori assumption on the stress tensor[END_REF][START_REF] Chapelle | The Finite Element Analysis of Shells -Fundamentals[END_REF] As a matter of fact, in shell models such Galerkin projections have been established to be asymptotically consistent with 3D elasticity, namely, the respective asymptotic analyses provide the same limit models. [START_REF] Chapelle | 3D-shell elements and their underlying mathematical model[END_REF] In addition, limit models obtained with polynomial assumptions are also interesting in relation to Galerkin discretizations of the variational problem, e.g. with finite elements used with a single element -of possibly high polynomial order -across the thickness.

If Assumption 4.2 holds, we can use Theorem 4.1 to show the existence and uniqueness of a solution of (4.4) because {u 0 + ξ

3 u 1 + (ξ 3 ) 2 u 2 , (u 0 , u 1 , u 2 ) ∈ V} is a closed subspace of V 3D . Theorem 4.2. For any ε = d L > 0, there exists a unique (u ε 0 , u ε 1 , u ε 2 ) ∈ V solution of (4.4).
We will now denote

u ε = u ε 0 + ξ 3 u ε 1 + (ξ 3 ) 2 u ε 2 ∈ V 3D
and study the behavior of this sequence when the parameter ε tends to zero.

Asymptotic analysis

In this section, we want to:

• identify a limit problem;

• show the existence and uniqueness of a solution to the limit problem;

• prove a preliminary weak convergence result;

• establish the strong convergence.

asymptoticAnaly-sis˙m3as

A surface-based electrophysiology model relying on asymptotic analysis 11

Limit problem

We start by considering the left-hand side of (4.4), i.e.

Ω σ • ∇u • ∇v dV = σ t Ω g • ∇u • ∇v dV + (σ l -σ t ) Ω τ ⊗ τ • ∇u • ∇v dV, = σ t Ω ∇u • ∇v dV + (σ l -σ t ) Ω τ • ∇u τ • ∇v dV, (4.5)
where we denote

u = u 0 + ξ 3 u 1 + (ξ 3 ) 2 u 2 and v = v 0 + ξ 3 v 1 + (ξ 3 ) 2 v 2 ,
and we will analyze the two integrals separately.

According to Proposition 3.4, we have

σ t ∇u• ∇v = σ t 2 α,β=1 ∂ α u(a αβ +ξ 3 ḡαβ )∂ β v +u 1 v 1 +2ξ 3 u 1 v 2 +u 2 v 1 +4 ξ 3 2 u 2 v 2 .
We can integrate all the polynomial terms in ξ 3 between -d 2 and + d 2 , and after integration the first term in the asymptotic expansion is given by

d σ t ω ∇u 0 • ∇v 0 + u 1 v 1 dS.
We can then denote by A l the bilinear form

A l (u 0 , u 1 ), (v 0 , v 1 ) = σ t ω ∇u 0 • ∇v 0 + u 1 v 1 dS, ∀(u 0 , u 1 ), (v 0 , v 1 ) ∈ H 1 (S) × L 2 (S).
Concerning next the second term in (4.5) we have

τ • ∇u τ • ∇v = 3 i,j=1 ( τ • g i )( τ • g j )∂ i u∂ j v = 2 α,β=1 ( τ • g α )( τ • g β )∂ α u∂ β v,
according to Assumption 4.2. Using the geometric definitions, we have that

g α (ξ 1 , ξ 2 , ξ 3 ) = a α (ξ 1 , ξ 2 ) + ξ 3 H α (ξ 1 , ξ 2 , ξ 3 ) with H α ∈ L ∞ (B), so that we can decompose ( τ • g α )( τ • g β ) into ( τ • g α )( τ • g β ) = (τ 0 • a α )(τ 0 • a β ) + ξ 3 J (ξ 1 , ξ 2 , ξ 3 ) cos 2 2θξ 3 d + (τ ⊥ 0 • a α )(τ ⊥ 0 • a β ) + ξ 3 J ⊥ (ξ 1 , ξ 2 , ξ 3 ) sin 2 2θξ 3 d + (τ 0 • a α )(τ ⊥ 0 • a β ) + (τ ⊥ 0 • a α )(τ 0 • a β ) +ξ 3 J (ξ 1 , ξ 2 , ξ 3 ) sin 2θξ 3 d cos 2θξ 3 d ,
where J , J ⊥ et J are all L ∞ (B). When integrating over ξ 3 , the following integrals appear

I n (θ) = 1 d n+1 d 2 -d 2 (ξ 3 ) n cos 2 2θξ 3 d dξ 3 = 1 2 n+1 1 -1 t n cos 2 (tθ)dt, J n (θ) = 1 d n+1 d 2 -d 2 (ξ 3 ) n sin 2 2θξ 3 d dξ 3 = 1 2 n+1 1 -1 t n sin 2 (tθ)dt, L n (θ) = 1 d n+1 d 2 -d 2 (ξ 3 ) n cos 2θξ 3 d sin 2θξ 3 d dξ 3 = 1 2 n+1 1 -1
t n cos(tθ) sin(tθ)dt.

We can now as before identify the first term of the asymptotic expansion, which gives d A f l with

A f l (u 0 , v 0 ) = (σ l -σ t ) ω I 0 (θ) τ 0 ⊗ τ 0 + J 0 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 •∇u 0 • ∇v 0 dS, where I 0 (θ) = 1 2 + 1 4θ sin(2θ) and J 0 (θ) = 1 2 - 1 4θ
sin(2θ) = 1 -I 0 (θ). Note that these functions can be extended by continuity in 0, with values 1 and 0, respectively. Remark 4.7.

• I 0 (θ) > 0, J 0 (θ) ≥ 0, ∀θ ∈ R, • J 0 (θ) = 0 (and I 0 (θ) = 1) if and only if θ = 0. In this case τ (ξ 1 , ξ 2 , ξ 3 ) = τ 0 (ξ 1 , ξ 2 ), i.e. the fiber direction is independent of the thickness.

Using similar manipulations, the first term in the asymptotic expansion of the right-hand side of (4.4) is given by dF 0 , with

F 0 (v 0 ) = ω f 0 v 0 dS, ∀v 0 ∈ H 1 (S). Defining V lim = H 1 (S) ∩ (BC ) × L 2 (S), our candidate limit problem finally consists in seeking (u l 0 , u l 1 ) ∈ V lim such that A l (u l 0 , u l 1 ), (v 0 , v 1 ) + A f l (u l 0 , v 0 ) = F 0 (v 0 ), ∀(v 0 , v 1 ) ∈ V lim . (4.6)
Theorem 4.3. There exists a unique (u l 0 , u l 1 ) ∈ V lim solution of (4.6).

Proof. It is straightforward to see that the bilinear form in the left-hand side of (4.6) is continuous and coercive on V lim , and symmetric, hence we can directly apply the Lax-Milgram theorem.

A surface-based electrophysiology model relying on asymptotic analysis 13

Note that, choosing v 0 = 0 in (4.6), we can show that u l 1 = 0, and

u l 0 ∈ H 1 (S) ∩ (BC ) satisfies, ∀v 0 ∈ H 1 (S) ∩ (BC ), σ t ω ∇u l 0 • ∇v 0 dS + (σ l -σ t ) ω I 0 (θ) τ 0 ⊗ τ 0 + J 0 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 •∇u l 0 • ∇v 0 dS = ω f 0 v 0 dS. (4.7)

Weak convergence

Theorem 4.4. We have the following convergences:

• (u ε 0 , u ε 1 ) converges weakly to (u l 0 , 0) in V lim when ε tends to 0, • εu ε 1 and ε 2 u ε 2 converge weakly to 0 in H 1 (S) when ε tends to 0.

The proof of Theorem 4.4 is divided into two steps:

(i) The first step consists in proving that a subsequence of (u ε 0 , u ε 1 ) ε converges weakly to a limit denoted by (u m 0 , u m 1 ) ∈ V lim when ε tends to 0. (ii) In the second step, we show that u m 0 = u l 0 , u m 1 = 0 and we can infer the theorem.

(i) Weak convergence of a subsequence Lemma 4.1. There exists one strictly positive constant

C independent of ε such that ∀(v 0 , v 1 , v 2 ) ∈ V, A 3D (v, v) ≥ C ε ∇v 0 L 2 (S) + v 1 L 2 (S) 2 + ε 3 ∇v 1 L 2 (S) + v 2 L 2 (S) 2 + ε 5 ∇v 2 2 L 2 (S) .
Lemma 4.2. There exists one strictly positive constant

C independent of ε such that ∀(v 0 , v 1 , v 2 ) ∈ V, F 3D (v) ≤ C ε ∇v 0 L 2 (S) + v 1 L 2 (S) + ε 2 ∇v 1 L 2 (S) + v 2 L 2 (S) + ε 3 ∇v 2 L 2 (S) .
We directly have Lemma 4.2. We must prove Lemma 4.1.

Proof. Using Propositions 3.2 and 3.3, we have

∀(v 0 , v 1 , v 2 ) ∈ V, A 3D (v, v) = σ t Ω ∇v • ∇v dV + (σ l -σ t ) Ω τ ⊗ τ • ∇v • ∇v dV ≥ σ t Ω ∇v • ∇v dV ≥ γ σ t Ω ∇v • ∇v √ a dξ ≥ γ σ t Ω αβ=1,2 g αβ ∂ α v∂ β v + v 2 1 + 4ξ 3 v 1 v 2 + 4 ξ 3 2 v 2 2 √ a dξ ≥ c Ω αβ=1,2 a αβ ∂ α v∂ β v + v 2 1 + 4ξ 3 v 1 v 2 + 4 ξ 3 2 v 2 2 √ a dξ ≥ c Ω ∇v • ∇v + v 2 1 + 4ξ 3 v 1 v 2 + 4 ξ 3 2 v 2 2 √ a dξ,
with Ω dependent on ε. We use the decomposition of v and we integrate over ξ 3

A 3D (v, v) ≥ c ω ε(∇v 0 • ∇v 0 + v 2 1 ) + ε 3 12 (∇v 1 • ∇v 1 + 4v 2 2 ) + ε 3 6 ∇v 0 • ∇v 2 + ε 5 80 ∇v 2 • ∇v 2 dS.
Using the Young inequality on the cross term, we have for almost all (ξ 1 , ξ 2 ) ∈ ω, and for all r > 0,

2∇v 0 (ξ 1 , ξ 2 ) • ∇v 2 (ξ 1 , ξ 2 ) ≥ - 1 r ∇v 0 (ξ 1 , ξ 2 ) • ∇v 0 (ξ 1 , ξ 2 ) -r∇v 2 (ξ 1 , ξ 2 ) • ∇v 2 (ξ 1 , ξ 2 ). With r = 1 10 ε 2 , we obtain ∀(v 0 , v 1 , v 2 ) ∈ V, A 3D (v, v) ≥ c ω ε 6 ∇v 0 • ∇v 0 + εv 2 1 + ε 3 12 (∇v 1 • ∇v 1 + 4v 2 2 ) + ε 5 240 ∇v 2 • ∇v 2 dS ≥ C ε ∇v 0 L 2 (S) + v 1 L 2 (S) 2 + ε 3 ∇v 1 L 2 (S) + v 2 L 2 (S) 2 +ε 5 ∇v 2 2 L 2 (S) .
These two lemmas directly allow to show the following result.

Lemma 4.3.

There exists a strictly positive constant C independent of ε such that

∇u ε 0 L 2 (S) + u ε 1 L 2 (S) + ε ∇u ε 1 L 2 (S) + u ε 2 L 2 (S) + ε 2 ∇u ε 2 L 2 (S) ≤ C.
Using Lemma 4.3, we can quickly finish the first step. According to the Poincaré inequality (Γ 1 of non-zero measure) and Lemma 4.3, u ε 0 , εu ε 1 et ε 2 u ε 2 are uniformly bounded in the H 1 (S)-norm and u ε 1 and εu ε 2 are uniformly bounded in the L 2 (S)norm. We can then extract a subsequence (also denoted by ε) such that there exists (u m 0 , u m 1 , ũm

1 , u m 2 , ũm 2 ) ∈ H 1 (S) × L 2 (S) × H 1 (S) × L 2 (S) × H 1 (S) such that

asymptoticAnaly-sis˙m3as

A surface-based electrophysiology model relying on asymptotic analysis 15

(u ε 0 , u ε 1 , εu ε 1 , εu ε 2 , ε 2 u ε 2 ) converges weakly to (u m 0 , u m 1 , ũm 1 , u m 2 , ũm 2 ) in H 1 (S)×L 2 (S)× H 1 (S) × L 2 (S) × H 1 (S).
Remark 4.8. We can show that ũm 1 = 0, since by compact injection εu ε

1 L 2 (S)
→ ũm 1 . Likewise, ũm 2 = 0.

(ii) Convergence to the solution of the limit problem Lemma 4.4.

∀(v 0 , v 1 , v 2 ) ∈ V, • 1 d A 3D (u ε , v) = A l (u ε 0 , u ε 1 ), (v 0 , v 1 ) + A f l (u ε 0 , v 0 ) + O(ε), • 1 d F 3D (v) = F 0 (v 0 ) + O(ε).
With Lemma 4.4 we can quickly conclude. According to Theorem 4.2, indeed, we have ∀v ∈ V, A 3D (u ε , v) = F 3D (v) and Lemma 4.4 gives

A l (u ε 0 , u ε 1 ), (v 0 , v 1 ) + A f l (u ε 0 , v 0 ) = F 0 (v 0 ) + O(ε), ∀(v 0 , v 1 ) ∈ V lim . (4.8)
Making ε tend to 0 in (4.8) and according to the first step, we have

A l (u m 0 , u m 1 ), (v 0 , v 1 ) + A f l (u m 0 , v 0 ) = F 0 (v 0 ), ∀(v 0 , v 1 ) ∈ V lim .
Using Theorem 4.3, we obtain that (u m 0 , u m 1 ) = (u l 0 , 0) where u l 0 verifies Problem (4.7). The limits (u m 0 , u m 1 , ũm 1 , ũm 2 ) = (u l 0 , 0, 0, 0) are independent of the subsequence so that we can infer the theorem. To complete the proof of Theorem 4.4, we must now show Lemma 4.4.

Proof.

A 3D (u ε , v) = σ t Ω ∇u ε • ∇v dV +(σ l -σ t ) Ω τ ⊗ τ • ∇u ε • ∇v dV, ∀(v 0 , v 1 , v 2 ) ∈ V,
with Ω dependent on ε. We start by studying the first term of A 3D (u ε , v), i.e.

σ t Ω ∇u ε • ∇v dV = σ t Ω ∇u ε • ∇v √ a dξ + σ t Ω ξ 3 ∇u ε • ∇v √ a(-2H + Kξ 3 )dξ.
If we can show that σ t Ω ∇u ε • ∇v √ a dξ converges when ε tends to 0, then

σ t Ω ξ 3 ∇u ε • ∇v √ a(-2H + Kξ 3
)dξ converges to 0 when ε tends to 0, so that we just need to study the first term.

σ t Ω ∇u ε • ∇v √ a dξ = σ t Ω αβ=1,2 g αβ ∂ α u ε ∂ β v + u ε 1 v 1 + 2ξ 3 u ε v 2 + u ε 2 v 1 + 4 ξ 3 2 u ε 2 v 2 √ a dξ = σ t Ω αβ=1,2 a αβ ∂ α u ε ∂ β v + u ε 1 v 1 + 2ξ 3 u ε 1 v 2 + u ε 2 v 1 + 4 ξ 3 2 u ε 2 v 2 √ a dξ + σ t Ω αβ=1,2 ξ 3 ḡαβ ∂ α u ε ∂ β v √ a dξ.
For the same reason, we just need to show that σ t

Ω ∇u ε • ∇v + u ε 1 v 1 + 2ξ 3 u ε 1 v 2 + u ε 2 v 1 + 4 ξ 3 2 u ε 2 v 2 √
a dξ converges when ε tends to 0.

We can apply the same method with the second term of A 3D (u ε , v) and we infer that we only need to show that (σ lσ t )

Ω τ ⊗ τ • ∇u ε • ∇v √ a dξ converges when ε tends to 0. We define A 3D 0 by ∀(u 0 , u 1 , u 2 ), (v 0 , v 1 , v 2 ) ∈ V, A 3D 0 (u, v) = σ t Ω ∇u • ∇v + u 1 v 1 + 2ξ 3 u 1 v 2 + u 2 v 1 + 4 ξ 3 2 u 2 v 2 √ a dξ + (σ l -σ t ) Ω τ ⊗ τ • ∇u • ∇v √ a dξ 3 ,
and we want to show that

1 d A 3D 0 (u ε , v) = A l ((u ε 0 , u ε 1 ), (v 0 , v 1 )) + A f l (u ε 0 , v 0 ) + O(ε), ∀(v 0 , v 1 , v 2 ) ∈ V.
We decompose

A 3D 0 (u ε , v) as A 3D 0 (u ε , v) = 4 n=0 I n (u ε , v) with 1 ε I 0 (u ε , v) = L A l (u ε 0 , u ε 1 ), (v 0 , v 1 ) + A f l (u ε 0 , v 0 ) , 1 ε I 1 (u ε , v) = εL 2 (σ l -σ t ) ω L 1 (θ) τ ⊥ 0 ⊗ τ 0 • ∇u ε 0 • ∇v 1 + τ 0 ⊗ τ ⊥ 0 • ∇u ε 0 • ∇v 1 dS + εL 2 (σ l -σ t ) ω L 1 (θ) τ ⊥ 0 ⊗ τ 0 • ∇u ε 1 • ∇v 0 + τ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇v 0 dS, 1 ε I 2 (u ε , v) = 1 12 ε 2 L 3 σ t ω ∇u ε 0 • ∇v 2 + ∇u ε 2 • ∇v 0 + ∇u ε 1 • ∇v 1 + 4u ε 2 v 2 dS + ε 2 L 3 (σ l -σ t ) ω I 2 (θ) τ 0 ⊗ τ 0 • ∇u ε 0 • ∇v 2 + J 2 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 0 • ∇v 2 dS + ε 2 L 3 (σ l -σ t ) ω I 2 (θ) τ 0 ⊗ τ 0 • ∇u ε 2 • ∇v 0 + J 2 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇v 0 dS + ε 2 L 3 (σ l -σ t ) ω I 2 (θ) τ 0 ⊗ τ 0 • ∇u ε 1 • ∇v 1 + J 2 (θ) τ ⊥ 0 ⊗ ⊥ 0 • ∇u ε 1 • ∇v 1 dS, 1 ε I 3 (u ε , v) = ε 3 L 4 (σ l -σ t ) ω L 3 (θ) τ ⊥ 0 ⊗ τ 0 • ∇u ε 1 • ∇v 2 + τ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇v 2 dS + ε 3 L 4 (σ l -σ t ) ω L 3 (θ) τ ⊥ 0 ⊗ τ 0 • ∇u ε 2 • ∇v 1 + τ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇v 1 dS, 1 ε I 4 (u ε , v) = 1 80 ε 4 L 5 (σ l -σ t ) ω ∇u ε 2 • ∇v 2 dS + ε 4 L 5 (σ l -σ t ) ω I 4 (θ) τ 0 ⊗ τ 0 • ∇u ε 2 • ∇v 2 + J 4 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇v 2 dS.
We have that

• ε ω L 1 (θ) τ ⊥ 0 ⊗ τ 0 • ∇u ε 0 • ∇v 1 + τ 0 ⊗ τ ⊥ 0 • ∇u ε 0 • ∇v 1 dS → ε→0 0, because u ε 0 H 1 (S) is bounded, • ε ω L 1 (θ) τ ⊥ 0 ⊗τ 0 •∇u ε 1 •∇v 0 + τ 0 ⊗τ ⊥ 0 •∇u ε 1 •∇v 0 dS → ε→0 0, because ũm 1 = 0.
We infer that

1 ε I 1 (u ε , v) → ε→0 0, ∀(v 0 , v 1 , v 2 ) ∈ V. With similar arguments, we can show that 1 ε I 2 (u ε , v), 1 ε I 3 (u ε , v) and 1 ε I 4 (u ε , v) all converge to 0 when ε tends to 0, for all (v 0 , v 1 , v 2 ) ∈ V. We clearly have 1 ε F 3D (v) → ε→0 F 0 (v 0 ), ∀(v 0 , v 1 , v 2 ) ∈ V, and
this concludes the proof of the theorem.

Strong convergence

We will complete our asymptotic analysis by proving that the convergences established in the previous section also hold -in fact -in the strong sense. Defining the two bilinear forms

Ãf l (ρ, µ) = (σ l -σ t ) ω I 0 (θ) τ 0 ⊗ τ 0 • ρ • µ dS, Āf l (ρ, µ) = (σ l -σ t ) ω J 0 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ρ • µ dS,
and the rotation operator

α : Vect( a 1 , a 2 ) -→ Vect( a 1 , a 2 ) τ -→ a 3 × τ
we will show that the following sum of perfect squares

D ε = A l u ε 0 + 1 12 d 2 u ε 2 -u l 0 , u ε 1 2 + Ãf l ∇ u ε 0 + I 2 I 0 d 2 u ε 2 - L 1 I 0 d α ∇u ε 1 -∇u l 0 2 + Āf l ∇ u ε 0 + I 2 J 0 d 2 u ε 2 + L 1 J 0 d α ∇u ε 1 -∇u l 0 2
, tends to zero with ε, where we denote A l (u, v) 2 = A l (u, v), (u, v) , and similarly for Ãf l and Āf l . In fact, D ε is constructed so that the expression

D ε = A l u ε 0 + 1 12 d 2 u ε 2 , u ε 1 2 + Ãf l ∇ u ε 0 + I 2 I 0 d 2 u ε 2 - L 1 I 0 d α ∇u ε 1 2 + Āf l ∇ u ε 0 + I 2 J 0 d 2 u ε 2 + L 1 J 0 d α ∇u ε 1 2
, appearing when developing D ε , gathers the main lower-order terms in the expansion A 3D (u ε , u ε ) as expressed in the following lemma, shown like in the proof of Lemma 4.4.

Lemma 4.5.

A 3D (u ε , u ε ) = dσ t ω ∇u ε 0 • ∇u ε 0 + (u ε 1 ) 2 dS + d 2 12 ω 2∇u ε 0 • ∇u ε 2 + ∇u ε 1 • ∇u ε 1 + (u ε 2 ) 2 dS + d 4 80 ω ∇u ε 2 • ∇u ε 2 dS + d(σ l -σ t ) ω I 0 (θ) τ 0 ⊗ τ 0 • ∇u ε 0 • ∇u ε 0 + J 0 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 0 • ∇u ε 0 dS + d(σ l -σ t ) d ω 2L 1 (θ) τ 0 ⊗ τ ⊥ 0 • ∇u ε 0 • ∇u ε 1 + 2L 1 (θ) τ ⊥ 0 ⊗ τ 0 • ∇u ε 0 • ∇u ε 1 dS + d(σ l -σ t ) d 2 ω 2I 2 (θ) τ 0 ⊗ τ 0 • ∇u ε 0 • ∇u ε 2 + 2J 2 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 0 • ∇u ε 2 dS + d(σ l -σ t ) d 2 ω I 2 (θ) τ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 1 + J 2 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 1 dS + d(σ l -σ t )L d 3 ω 2L 3 (θ) τ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 2 + 2L 3 (θ) τ ⊥ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 2 dS + d(σ l -σ t )L d 4 ω I 4 (θ) τ 0 ⊗ τ 0 • ∇u ε 2 • ∇u ε 2 + J 4 (θ) τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇u ε 2 dS + R(u ε ),
where the degree of ξ 3 in the terms of We can now state and prove our final convergence result. Theorem 4.5. We have the following convergences:

1 d R(u ε ) = 1 Lε R(u ε
• (u ε 0 , u ε 1 )
converges strongly to (u l 0 , 0) in V lim when ε tends to 0, • εu ε 1 and ε 2 u ε 2 converge strongly to 0 in H 1 (S) when ε tends to 0.

In addition, εu ε 2 converges strongly to 0 in L 2 (S) when ε tends to 0.

Proof. Defining Dε = 1 12 d 2 σ t ω ∇u ε 1 • ∇u ε 1 + (u ε 2 ) 2 dS + 1 180 d 4 σ t ω ∇u ε 2 • ∇u ε 2 dS + d 2 (σ l -σ t ) ω I 2 - L 2 1 J 0 τ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 1 dS + ω J 2 - L 2 1 I 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 1 dS + 2d 3 (σ l -σ t ) ω L 3 -L 1 I 2 I 0 τ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 2 dS + ω L 3 -L 1 J 2 J 0 τ ⊥ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 2 dS + d 4 (σ l -σ t ) ω I 4 - I 2 2 I 0 τ 0 ⊗ τ 0 • ∇u ε 2 • ∇u ε 2 dS + ω J 4 - J 2 2 J 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇u ε 2 dS ,
we will decompose the proof into two steps:

• in the first step, we show that D ε + Dε → ε→0 0,

• in the second step, we show that Dε ≥ η εu ε 1 2

H 1 (S) + εu ε 2 2 L 2 (S) + ε 2 u ε 2 2 H 1 (S) , with η > 0.
The conclusion is then at hand, because it implies that D ε → ε→0 0 and Dε → ε→0 0, and in particular we have the following convergences

A l u ε 0 + 1 3 d 2 u ε 2 -u l 0 , u ε 1 2 → ε→0 0, εu ε 1 H 1 (S) → ε→0 0, εu ε 2 L 2 (S) → ε→0 0, ε 2 u ε 2 H 1 (S) → ε→0 0.
We infer

u ε 0 H 1 (S) -→ ε→0 u l 0 , u ε 1 L 2 (S) -→ ε→0 0, εu ε 1 H 1 (S) -→ ε→0 0, εu ε 2 L 2 (S) -→ ε→0 0, ε 2 u ε 2 H 1 (S) -→ ε→0 0.
First step We can develop D ε as

D ε = D ε -2A l u ε 0 + 1 12 d 2 u ε 2 , 0 , (u l 0 , 0) -2 Ãf l ∇ u ε 0 + I 2 I 0 d 2 u ε 2 - L 1 I 0 d α ∇u ε 1 , ∇u l 0 -2 Āf l ∇ u ε 0 + I 2 J 0 d 2 u ε 2 + L 1 J 0 d α ∇u ε 1 , ∇u l 0 + A l (u l 0 , 0), (u l 0 , 0) + Ãf l ∇u l 0 , ∇u l 0 + Āf l ∇u l 0 , ∇u l 0 , We remark that Ãf l ∇u l 0 , ∇u l 0 + Āf l ∇u l 0 , ∇u l 0 = A f l u l 0 , u l 0 ,
and according to Theorem 4.3, we have A l (u l 0 , 0), (u l 0 , 0) + A f l u l 0 , u l 0 = F 0 (u l 0 ). Furthermore, according to Theorem 4.4, we have

A l u ε 0 + 1 3 d 2 u ε 2 , 0 , (u l 0 , 0) + Ãf l ∇ u ε 0 + I 2 I 0 d 2 u ε 2 - L 1 I 0 d α ∇u ε 1 , ∇u l 0 + Āf l ∇ u ε 0 + I 2 J 0 d 2 u ε 2 + L 1 J 0 d α ∇u ε 1 , ∇u l 0 → ε→0 A l u l 0 , u l 0 + A f l u l 0 , u l 0 = F 0 (u l 0 ). We infer D ε = D ε -F 0 (u l 0 ) + O(ε). (4.9)
Developing D ε and using Lemma 4.5, we have

D ε = 1 d A 3D (u ε , u ε ) - 1 d R(u ε ) - 1 12 d 2 σ t ω ∇u ε 1 • ∇u ε 1 + (u ε 2 ) 2 dS - 1 180 d 4 σ t ω ∇u ε 2 • ∇u ε 2 dS -d 2 (σ l -σ t ) ω I 2 - L 2 1 J 0 τ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 1 dS + ω J 2 - L 2 1 I 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 1 dS -2d 3 (σ l -σ t ) ω L 3 -L 1 I 2 I 0 τ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 2 dS + ω L 3 -L 1 J 2 J 0 τ ⊥ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 2 dS -d 4 (σ l -σ t ) ω I 4 - I 2 2 I 0 τ 0 ⊗ τ 0 • ∇u ε 2 • ∇u ε 2 dS + ω J 4 - J 2 2 J 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇u ε 2 dS = 1 d F 3D (u ε ) - 1 d R(u ε ) -Dε = F 0 (u l 0 ) -Dε + O(ε).
The equation (4.9) then becomes D ε + Dε = O(ε), which concludes the first step.

H 1 (S) + εu ε 2 2 L 2 (S) + ε 2 u ε 2 2 H 1 (S) . Dε = 1 12 d 2 σ t ω ∇u ε 1 • ∇u ε 1 + (u ε 2 ) 2 dS + 1 180 d 4 σ t ω ∇u ε 2 • ∇u ε 2 dS +d 2 (σ l -σ t ) ω I 2 - L 2 1 J 0 τ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 1 dS + ω J 2 - L 2 1 I 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 1 dS +2d 3 (σ l -σ t ) ω L 3 -L 1 I 2 I 0 τ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 2 dS + ω L 3 -L 1 J 2 J 0 τ ⊥ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 2 dS +d 4 (σ l -σ t ) ω I 4 - I 2 2 I 0 τ 0 ⊗ τ 0 • ∇u ε 2 • ∇u ε 2 dS + ω J 4 - J 2 2 J 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇u ε 2 dS .
We have

1 3 d 2 σ t ω ∇u ε 1 • ∇u ε 1 + (u ε 2 ) 2 dS + 4 45 d 4 σ t ω ∇u ε 2 • ∇u ε 2 dS ≥ η εu ε 1 2 H 1 (S) + εu ε 2 2 L 2 (S) + ε 2 u ε 2 2 H 1 (S) ,
so that the result holds if we can show the positiveness of Dε defined by

Dε = (σ l -σ t ) d 2 ω I 2 - L 2 1 J 0 τ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 1 dS +2d 3 ω L 3 -L 1 J 2 J 0 τ ⊥ 0 ⊗ τ 0 • ∇u ε 1 • ∇u ε 2 dS +(σ l -σ t )d 4 ω J 4 - J 2 2 J 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 2 • ∇u ε 2 dS +(σ l -σ t ) d 2 ω J 2 - L 2 1 I 0 τ ⊥ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 1 dS +2d 3 ω L 3 -L 1 I 2 I 0 τ 0 ⊗ τ ⊥ 0 • ∇u ε 1 • ∇u ε 2 dS +(σ l -σ t )d 4 ω I 4 - I 2 2 I 0 τ 0 ⊗ τ 0 • ∇u ε 2 • ∇u ε 2 dS.
We can decompose Dε into two sums of the form au 2 + 2cuv + bv 2 . We then use the property that au 2 + 2cuv + bv 2 is positive for all (u, v) if and only if a, b ≥ 0 and c 2 ≤ ab. Therefore, Dε is positive if and only if

I 2 - L 2 1 J 0 ≥ 0, J 4 - J 2 2 J 0 ≥ 0, L 3 -L 1 J 2 J 0 2 ≤ I 2 - L 2 1 J 0 . J 4 - J 2 2 J 0 , J 2 - L 2 1 I 0 ≥ 0, I 4 - I 2 2 I 0 ≥ 0, L 3 -L 1 I 2 I 0 2 ≤ J 2 - L 2 1 I 0 . I 4 - I 2 2 I 0 .
The inequalities in the left and center columns are simple consequences of the Cauchy-Schwarz inequality, while those in the right column can easily be checked analytically using, e.g., symbolic computation software.

Remark 4.9. Various possible extensions can be considered for the asymptotic setting assumptions:

• There is no particular difficulty in considering non-homogeneous Neumann boundary conditions on the lateral surface Γ 2 in the asymptotic analysis, with similar assumptions as for f regarding ξ 3 -regularity. • We conjecture that polynomial assumptions of higher degree could be handled in a similar manner and would provide similar convergence results to the same limits. • We could also consider more general forms of the angular variations of the fibers across the thickness of the type τ (ξ 1 , ξ 2 , ξ 3 ) = τ 0 (ξ 1 , ξ 2 ) cos θ ξ 1 , ξ 2 , 2ξ 3 /d + τ ⊥ 0 (ξ 1 , ξ 2 ) sin θ ξ 1 , ξ 2 , 2ξ 3 /d , which would lead to extended definitions of the geometric coefficients I i , J i and L i . The convergence proofs could then be performed verbatim up to the final chains of inequalities to be satisfied by these coefficients, which would have to be checked on a case-by-case basis, except for those directly following from the Cauchy-Schwarz inequality. Other types of scaling of the transverse coordinate could also be investigated in the definition of the fiber direction, e.g. with d α substituted for d, but our objective here was to obtain a finite variation of the angle across the thickness in view of the applications considered.

Surface-based bidomain model and numerical assessments

Following our above detailed asymptotic analysis, we directly adapt the result to propose the surface-based bidomain model, for all t > 0, find V m (•, t) ∈ H 1 (S), u e (•, t) ∈ H 1 (S) and w(•, t) ∈ L ∞ (S) with ω u e dS = 0, such that

                 A m ω C m ∂V m ∂t + I ion (V m , w) φ dS + ω σ i • ∇V m + ∇u e • ∇φ dS = A m ω I app (V m , w)φ dS, ω (σ i + σ e ) • ∇u e • ∇ψ dS + ω σ i • ∇V m • ∇ψ dS = 0, ∂ t w + g(V m , w) = 0, in ω × (0, T ).
for all φ, ψ ∈ H 1 (S) such that ω ψ dS = 0 and with σ i,e = σ t i,e a + (σ l i,eσ t i,e )(I 0 (θ)τ 0 ⊗ τ 0 + J 0 (θ)τ ⊥ 0 ⊗ τ ⊥ 0 ), where a denotes the surface-based metric tensor -namely, the restriction of the 3D metric tensor to vectors lying in the tangential plane -and recalling that τ 0 is a unit vector parallel to the local fiber direction on the midsurface, 2θ is the angle between the fiber directions on the lower and upper boundary surfaces -referred to as the endocardium and epicardium, respectively, in the heart -and I 0 (θ) = [START_REF] Ascher | Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations[END_REF] 2 + 1 4θ sin(2θ), J 0 (θ) = 1 2 -1 4θ sin(2θ). We will now perform a numerical assessment of our proposed model by comparing the results given by the 3D model (2.1) and the surface model (5.1) for two different geometries. For these comparisons, we limit ourselves to the Mitchell-Schaeffer model to avoid undue technicalities in the ionic model calibration and solution, although of course other ionic models could be considered, e.g. more adapted to the atria cells. The values of the parameters used in the simulations are given in Table 1. The tests are performed with the finite element library FELiScE, developed at Inria. We use P 1 -Lagrange finite elements, and as time discretization scheme a standard Backward Differentiation Formula (BDF) of order two, see e.g. Ref. 1. Although we primarily aim at a numerical assessment of the reduced model, we will use realistic values characteristic of atrial electrophysiology for all dimensions and parameters.

Planar test case

In this first test case, the 3D domain is a rectangular cuboid with dimensions 10 × 10 × 0.2 (all dimensions given in cm), hence the midsurface is a 10 × 10 square, which we discretize into 200 elements in each direction, namely, 40,401 vertices and 80,000 triangular elements. The 3D mesh is obtained by extrusion of the 2D mesh using 4 elements across the thickness, which gives 202,005 vertices and 960,000 tetrahedral elements. The conductivities are σ t i = 4.0 10 -4 , σ t e = 2.2 10 -3 , σ l i = 4.0 10 -3 , σ l e = 4.0 10 -3 , (all in S.cm -1 ). The fiber directions in the 3D mesh vary across the thickness only, and by an angle Therefore, on the midsurface the fiber direction is aligned with the vector

τ 0 = √ 2 2   1 1 0   ,
and we have θ = π 8 . The model is initialized at t = 0 with V m = -80 mV, u e = 0 mV, and we trigger the wave by applying a current I app = 0.5 mA.cm -2 in a circle (or cylindrical region in 3D) centered at (0 , 5) and of radius r = 0.55 during 1 ms.

In Figure 3, we compare the transmembrane potential given by the asymptotic surface model, the 3D model and a naive 2D model obtained by taking the value of the 3D conductivity tensor on the midsurface, namely, with σ i,e = σ t i,e a + (σ l i,eσ t i,e )τ 0 ⊗ τ 0 . We observe an excellent agreement between the asymptotic surface model and the 3D model, whereas the naive model exhibits a slower front wave in the direction perpendicular to the midsurface fiber. In addition, we plot in Figure 4 the difference |u ref,3Du 2D | on the midsurface. We notice that the differences between these two models are very limited and are narrowly concentrated near the wave front.

asymptoticAnaly-sis˙m3as

A surface-based electrophysiology model relying on asymptotic analysis 25 We also report in Table 2 the normalized l ∞ and l 2 differences between the asymptotic surface and the 3D models, i.e.

sup nodes |u ref,3D -u 2D | Vmax-Vmin and 1 Vmax-Vmin 1 N odes nodes |u ref,3D -u 2D | 2 1 2
, respectively. Note that V max -V min is a typical normalizing constant, very similar to the L ∞ -norm of u ref,3D and also meaningful from a modeling standpoint since it appears in the ionic model. The relative l ∞ errors appear to be quite high, but when closely studying Figures 3 and4 we see that this corresponds to very small shifts in the front location, which due to the dramatic steepness of the front induces significant -albeit very localized -errors. This is confirmed by the l 2 error values that are much smaller, indeed. We also point out that these errors are quite stable over time. Table 3 shows the same quantities when the asymptotic surface model is substituted by the naive 2D model. These results quantitatively confirm the superiority of the asymptotic model with respect to the naive model, and the excellent accuracy of the asymptotic model shows the relevance of the asymptotic analysis carried out in the previous sections when the fiber orientation rapidly varies across the thickness.

Spiral wave re-entry on a cylinder

This test case is motivated by so-called spiral waves, a fascinating phenomenon in cardiac electrophysiology, also very important in that it is often argued to be responsible for atrial or ventricular fibrillation, see e.g. Ref. 11. This complex behavior provides a challenging test case for our proposed model, and we will also consider a non-planar geometry, indeed.

The domain is a 3D half-cylinder, with internal radius r min = 4.9, external radius r max = 5.1 and height h = 10 (all dimensions again given in cm). Thus, the midsurface is a half-cylinder of radius 5.0 and height 10. The 3D mesh contains 105,080 vertices and 499,080 tetrahedral elements with 4 elements across the thickness, and the 2D mesh has 21,016 vertices and 41,590 triangular elements. The fiber directions in the 3D geometry vary across the thickness by an angle π 2 around the midsurface direction given by the vector

τ 0 = 1 5   -y x 0   ,
so that the fiber directions in the 2D mesh are characterized by this τ 0 and θ = π 4 . For this experiment, the conductivities are σ t i = 4.0 10 -4 , σ t e = 2.2 10 -3 , σ l i = 4.0 10 -3 , σ l e = 4.0 10 -3 (in S.cm -1 ). In order to initiate a spiral wave, we follow the procedure used in Ref. 9. The model is initialized at t = 0 with V m = -80 mV, u e = 0 mV, and we first generate a wave in the horizontal direction. Once the wave front has formed, it starts traveling along the cylinder (see Figure 5 at t = 80.0 ms). Just after the end of the repolarization (see Figure 5 at t = 300.0 ms), we depolarize a second region at time t = 375.0 ms. The new depolarization front revolves around the region that has just started repolarizing, which triggers a spiral wave.

We display in Figure 5 the results obtained with the 3D and surface models at successive times, and in Figure 6 plotted on the midsurface. As in the planar test case experiment, the differences between the 3D and asymptotic surface models are quite limited and very narrowly concentrated along the wave front. By contrast, the naive 2D model appears to be of poor accuracy, and in fact altogether destroys the spiral wave after 40 ms. This is fully confirmed by Tables 4 and5 showing the normalized differences in norms l ∞ and l 2 (less than 4% in l 2 -norm, here, for the asymptotic model). These results substantiate the validity of our proposed surface-based electrophysiology model to accurately represent such complex behaviors, and including over long time spans.

Fig. 1 .

 1 Fig. 1. Geometric domain and midsurface S descriptions

4 . 3 . 1 .

 431 Proposition Introducing the second and third fundamental forms respectively defined for all α, β = 1, 2 by

Remark 4 . 5 .

 45 Assumption 4.2 implies Assumption 4.1.

  ) is high enough to have the convergence to 0 when ε tends to 0. April 24, 2013 10:6 WSPC/INSTRUCTION FILE asymptoticAnaly-sis˙m3as A surface-based electrophysiology model relying on asymptotic analysis 19

A

  surface-based electrophysiology model relying on asymptotic analysis 21 Second step We want to show that Dε ≥ η εu ε

Fig. 3 .Fig. 4 .

 34 Fig. 3. Planar test case -Comparison of asymptotic surface model (left) 3D model (center) and 2D naive model (right) on the midsurface at 8 consecutive times

Fig. 5 .

 5 Fig. 5. Spiral wave on cylinder -Comparison of asymptotic surface model (left), 3D model (center) and naive 2D model (right) on the midsurface 8 consecutive times

Fig. 6 .

 6 Fig. 6. Spiral wave on cylinder -Point-wise difference between the 3D model and the asymptotic surface model on the midsurface (to be compared with the typical variation range of the transmembrane potential, i.e. about 100 mV)

Table 1 .

 1 Cell membrane parameters

	Am	Cm	τ in	τout	τopen	τ close	Vgate	V min	Vmax
	(cm -1 ) (mF.cm -2 ) (cm 2 .mA -1 ) (cm 2 .mA -1 )	(ms)	(ms)	(mV)	(mV)	(mV)
	200.0	10 -3	4.0	90.0	100.0	100.0	-67.0 -80.0	20.0

Table 2 .

 2 Planar test case -Error norms between the 3D and asymptotic surface models, on

	the midsurface
	T(ms)	30	50	60	70	100	140
	| • |∞	1.081 10 -1	1.582 10 -1	1.792 10 -1	2.364 10 -1	3.172 10 -1	3.687 10 -1
	| • | l 2	6.764 10 -3	1.179 10 -2	1.442 10 -2	1.624 10 -2	1.933 10 -2	1.934 10 -2

Table 3 .

 3 Planar test case -Error norms between the 3D and naive 2D models, on the midsurface

	T(ms)	30	50	60	70	100	140
	| • |∞	4.772 10 -1	7.447 10 -1	8.339 10 -1	9.281 10 -1	1.000	1.000
	| • | l 2	4.020 10	8.218 10 -2	1.035 10 -1	1.151 10 -1	1.530 10 -1	1.757 10 -1

Table 4 .

 4 Spiral wave on cylinder -Error norms between the 3D and asymptotic surface models, on the midsurface

	April 24, 2013 10:6 WSPC/INSTRUCTION FILE sis˙m3as	asymptoticAnaly-
		A surface-based electrophysiology model relying on asymptotic analysis 27
	T(ms)	80	300	390	700	800	900
	| • |∞	2.351 10 -1	1.467 10 -2	3.747 10 -2	1.525 10 -1	4.913 10 -1	4.229 10 -1
	| • | l 2	8.375 10 -3	2.315 10 -3	2.797 10 -3	2.019 10 -2	3.744 10 -2	3.228 10 -2

Table 5 .

 5 Spiral wave on cylinder -Error norms between the 3D and naive 2D models, on the midsurface

	T(ms)	80	300	390	700	800	900
	| • |∞	9.110 10 -1	1.026 10 -1	2.338 10 -1	9.548 10 -1	9.525 10 -1	8.788 10 -1
	| • | l 2	1.419 10 -1	4.645 10 -2	4.026 10 -2	4.169 10 -1	4.542 10 -1	4.318 10 -1

D. Chapelle, A. Collin and J.-F. Gerbeau

D. Chapelle, A. Collin and J.-F. Gerbeau

D. Chapelle, A. Collin and J.-F. Gerbeau

D. Chapelle, A. Collin and J.-F. Gerbeau

D. Chapelle, A. Collin and J.-F. Gerbeau

D. Chapelle, A. Collin and J.-F. Gerbeau