Yves Coudière

Jacques Henry

Simon Labarthe
email: simon.labarthe@u-bordeaux2.fr

A two layers monodomain model of cardiac electrophysiology of the atria : numerical illustration

Keywords: Cardiac modeling, Atrial model, Surface model, Asymptotic analysis. Subclass: MSC 92C30, MSC 92C50 and 35Q92. 1

. We investigate a two-dimensional model with two coupled, superimposed layers that allows to introduce three-dimensional heterogeneities, but retains a reasonable computational cost. We introduce the mathematical derivation of this model and error estimates with respect to the three-dimensional model. We give some numerical illustrations of its interest: we numerically show its convergence for vanishing thickness, introduce an optimization process of the coupling coefficient and assess its validity on physiologically relevant geometries. Our model would be an efficient tool to test the influence of three-dimensional fibre direction heterogeneities in reentries or atrial arrhythmia without using three-dimensional models.

Introduction

Modeling the electrophysiology of cardiac tissues is considered an investigation tool for clinical and fundamental research. Theoretical studies and numerical simulations have been recognized an efficient way to improve our knowledge of arrhythmia genesis [START_REF] Zemlin | Spontaneous onset of atrial fibrillation[END_REF][START_REF] Haissaguerre | Atrial fibrillatory cycle length: computer simulation and potential clinical importance[END_REF], perpetuation [START_REF] Johnw | Criterion for stable reentry in a ring of cardiac tissue[END_REF], and ablation efficiency [ROTTER et al., 2007[START_REF] Nagaiah | Optimal control approach to termination of re-entry waves in cardiac electrophysiology[END_REF] in atria.

In three-dimensional tissues, the propagation of the action potential is modeled by the monodomain or bidomain systems of equations. These are reactiondiffusion equations, where the reaction represents the ionic flow through the membrane of the cardiac cells that compose the tissue and the diffusion accounts for the diffusion of electrical charges.

Models of the propagation of the action potential through the atria are often formulated as monodomain or bidomain systems of equations on bi-dimensional manifolds [START_REF] Haissaguerre | Atrial fibrillatory cycle length: computer simulation and potential clinical importance[END_REF], ROTTER et al., 2007[START_REF] Dang | Evaluation of ablation patterns using a biophysical model of atrial fibrillation[END_REF]. Such surface models take advantage of the thinness of the atria with respect to the length scale of the heart, and drastically reduce the numerical cost of their resolution, hence allowing thorough in-silico investigations of atrial arrhythmias. These surface models approximate the first term of the asymptotic expansion (with respect to the thickness) of the three-dimensional equations, as presented in section 4. The thickness itself is not part of the final model, since it is assumed to vanish.

On the other hand, there exist structural and functional heterogeneities within the atria, such as fibre direction discontinuities through the wall. For instance, several superimposed layers with different fibre directions were described in the atrial wall or in the ostium zones of the pulmonary veins [START_REF] Ho | Architecture of the pulmonary veins: relevance to radiofrequency ablation[END_REF][START_REF] Sy Ho | Anatomy of the left atrium: implications for radiofrequency ablation of atrial fibrillation[END_REF][START_REF] Hocini | Electrical conduction in canine pulmonary veins: Electrophysiological and anatomic correlation[END_REF][START_REF] Saito | Left atrial myocardial extension onto pulmonary veins in humans: anatomic observations relevant for atrial arrhythmias[END_REF]. In these anatomical regions, the transition of fibres directions between layers can be very abrupt : consecutive fibre layers can have orthogonal directions.

These abrupt transitions can trigger complex propagation patterns [START_REF] Vetter | Epicardial fiber organization in swine right ventricle and its impact on propagation[END_REF] (and also Figure (3)), with several directions of anisotropy, and are suspected to be a substrate for re-entries or arrhythmia [START_REF] Hocini | Electrical conduction in canine pulmonary veins: Electrophysiological and anatomic correlation[END_REF], Ho et al., 2002[START_REF] Nattel | New ideas about atrial fibrillation 50 years on[END_REF].

Functional transmural heterogeneities have also been demonstrated, in particular during atrial arrhythmia [START_REF] Eckstein | Time course and mechanisms of endo-epicardial electrical dissociation during atrial fibrillation in the goat[END_REF]. Electrical dissociations have been observed in vivo during simultaneous endo and epicardial recording on dog heart [START_REF] Derakhchan | Method for simultaneous epicardial and endocardial mapping of in vivo canine heart: Application to atrial conduction properties and arrhythmia mechanisms[END_REF]. Differences in effective refractory period have also been measured on porcine [START_REF] Michowitz | Electrophysiological differences between the epicardium and the endocardium of the left atrium[END_REF] and goat [START_REF] Eckstein | Time course and mechanisms of endo-epicardial electrical dissociation during atrial fibrillation in the goat[END_REF] atrial fibrillation models.

The usual surface models do not account for such heterogeneities. The model proposed in [START_REF] Chapelle | A surfacebased electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling[END_REF] is a surface bidomain system of equations with anisotropic conductivity derived from fibre direction varying smoothly through the atrial wall. This derivation is based on a rigorous mathematical study, including curvature effects. Although it seems to be an improvement over the usual models, it still does not account for the complex propagation pattern mentioned above.

In order to overcome this difficulty, models with several surface layers have been proposed [START_REF] Weber | On the influence of a volume conductor on the orientation of currents in a thin cardiac issue[END_REF][START_REF] Jacquemet | A biophysical model of atrial fibrillation and electrograms : formulation, validation and applications[END_REF][START_REF] Gharaviri | A computer model of endo-epicardial electrical dissociation and transmural conduction during atrial fibrillation[END_REF]. The main idea is to include the structural heterogeneities of leading importance while keeping the convenience of a surface model. However, the mathematical foundations of these models have not been studied in depth, and the physiological assumptions on which they rely are still unclear.

In this paper, we show that a special attention must be paid to the physiological scaling of the equations. We notably show that the balance between the reaction terms and transverse diffusion (conductivity through the wall) entirely determines the regime of the solutions. If transverse diffusion dominates, then the usual surface models are sufficient to represent the average behavior of the action potential. But if reaction cannot be neglected compared to transverse diffusion, then complex patterns (Figures (2)-(3)) are observed, and the transmural effects must not be neglected. This is typically the case during electrical depolarization for human atria, as discussed in section 3.2. Although remaining in the framework of vanishing thickness, we try to account for this regime by introducing the second order term in the asymptotic expansion of the solutions, instead of simply the first one. But this would not yield a model that clearly differentiates between different layers. Hence we propose a model that represents the evolution of the transverse averages of the transmembrane voltage in each of the tissue's layers. Consequently, we design a surface model with two layers having distinct conductivity tensor and electrophysiology source function. With this new vision, we can see how the thickness appears in the equations, and we account for the main structural and functional observations described above.

For a slab of tissue with two distinct layers, the model is formally derived from the three-dimensional monodomain equations in two steps: the second order term of the expansion of the 3D model with respect to the thickness is first computed. Then the averages of this second order model are shown to solve a set of two coupled surface monodomain equations, associated to the two layers of the tissue. The error between this two-surfaces model and the transverse averages of the threedimensional model is of order ε 3 , where ε is the aspect ratio of the atria. This result is demonstrated in an other paper [START_REF] Coudière | An asymptotic two layers monodomain model of cardiac electrophysiology in the atria[END_REF]. This theoretical approach for vanishing thickness is supplemented by a numerical procedure that improves the two layer model for physiological thickness of tissue. The convergence of the two layer model is numerically illustrated and its validity is assessed on physiologically relevant geometries.

The paper is organized as follows. Section 2 briefly recalls the monodomain equations. Section 3 introduces the geometrical framework of our paper, defines the small parameter ε, explains the scaling of the equations and discusses the different regimes observed. It justifies the fact that the thickness cannot be neglected during the depolarization process in the atria. Section 4 presents a derivation of the second order surface model similar (but one order more accurate) to the usual surface models, and introduces convergence results in 6. Section 5 proposes a derivation of our two layers model with two systems of monodomain equations formulated on a surface, and introduces convergence results. Section 7 is devoted to the numerical illustrations: it numerically shows the convergence of the surface models, formalizes an optimization method of the coupling coefficient and address the robustness of this model for geometries with physiological characteristics. The conclusions are discussed in section 8.

The monodomain equations

The bidomain equations, first introduced in [START_REF] Clerc | Directional differences of impulse spread in trabecular muscle from mammalian heart[END_REF][START_REF] Tung | A Bi-Domain Model for describing Ischemic Myocardial D-C Potentials[END_REF], are a degenerate system of two anisotropic reaction-diffusion equations coupled to a set of ordinary differential equations. They describe the evolution of the intra-and extra-cellular potential at the tissue scale [START_REF] Krassowska | Homogenization of syncytial tissues[END_REF], and they can mimic complex phenomena such as virtual electrodes [START_REF] Sepulveda | Current injection into a twodimensional anisotropic bidomain[END_REF]. Under the equal anisotropy ratio assumption, the system can be reduced to the single monodomain reaction-diffusion equation. In the vast majority of applications, the solution to this monodomain equation is very close to the solution to the bidomain equation [START_REF] Potse | A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF]. The models proposed and discussed in this paper are based on the monodomain approximation although it is straightforward to derive a bidomain version of them.

We set the monodomain equations on an open subset Ω of R 3 and we denote by (u, w) its solution. The system reads

A (C∂ t u + f (u, w)) = div(σ ∇u) in (0, +∞) × Ω, (1)
∂ t w + g(u, w) = 0 in (0, +∞) × Ω, (2)
where A, C and σ are physiological parameters that refer to, respectively, the ratio of surface of membrane to total volume (in cm -1 ), the membrane capacitance (in µFcm -2 ) and the conductivity (in mScm -1 ). The unknowns u(t, x) ∈ R (in mV) describes the transmembrane potential and the m variables w(t, x) ∈ R m model the electrophysiological membrane dynamics. The functions g : R × R m → R m and f : R × R m → R define the dynamics of the conductance of the ionic channels through the membrane and of the ionic currents flowing through that channels.

The myocardial domain is represented by the domain Ω. In this domain, the cardiac tissue is composed of fibres, that in turn are structured into laminae. The tissue is modeled as a homogenized continuum with heterogeneous and anisotropic electrical conductivity σ :

∀x ∈ Ω, σ (x) = 3 ∑ i=1 σ i ν i (x)ν T i (x),
where (ν 1 , ν 2 , ν 3 ) is the orthonormal basis in R 3 aligned on the fibre and laminae directions and 0 < σ := σ 3 ≤ σ 2 ≤ σ 1 := σ are constant parameters. The diffusion operator is then bounded and uniformly elliptic:

∀x ∈ Ω, σ |ξ | 2 ≤ ξ T σ (x)ξ = 3 ∑ i=1 σ i |ξ i | 2 ≤ σ |ξ | 2 .
Equations ( 1) and ( 2) are supplemented with the Neumann boundary conditions

σ ∇u • n = 0 in (0, +∞) × ∂ Ω (3)
that models a complete electrical insulation of the heart. We add the initial conditions u(0, x) = u 0 (x), w(0, x) = w 0 (x) a.e. x ∈ Ω.

(4)

We assume that the monodomain problem (1) and (2) with the boundary and initial conditions (3) and ( 4) is well-posed. We then suppose that there is unique solutions u and w defined for a.e. x ∈ Ω and for all t > 0. Furthermore, we assume that this solution is regular.

3 Fibre distribution and dimensional analysis

A two-layers slab of myocardium

An ideal slab of cardiac tissue Ω is considered. It is composed by two superimposed layers with the same thickness h > 0 but different fibre directions. We set Ω = ω × (-h, h) ⊂ R 3 , where ω is an open bounded subset of R 2 , and the two layers Ω (1) = ω × (0, h) and Ω (2) = ω × (-h, 0). We decompose the boundary of each layer into the external boundary and the interface between the layers, namely

Γ (k) = ∂ Ω (k) ∩ ∂ Ω for k = 1, 2 and Σ = ω × {0}.
We express the coordinates of a point x ∈ Ω as x = (x , z) with x ∈ ω and -h < z < h. We assume that the fibre distribution is homogeneous through each layer. For k = 1, 2, let θ (k) be the angle of the fibres with the first coordinate axis. In each layer we build an orthonormal basis ν (k) i for i = 1, 2, 3 respectively along the fibre direction, orthogonal to the fibres in the lamina plane and orthogonal to the lamina by

ν (k) 1 =   cos(θ (k) ) sin(θ (k) ) 0   , ν (k) 2 =   -sin(θ (k) ) cos(θ (k) ) 0   , ν (k) 3 =   0 0 1   , (5)
where θ (1) and θ (2) are fixed angles. Then, the conductivity tensors in the layers read

σ (k) = 2 ∑ j=1 σ (k) j ν (k) j ν (k) j T =σ (k) +σ (k) 3 ν (k) 3 ν (k) 3 T for k = 1, 2. ( 6
)
where σ (k) is the two dimensional conductivity tensor of the layer number k.

The conductivity tensor σ in Ω is then defined by σ 2) . Hence, for initial data u 0 and w 0 in Ω, the monodomain equations ( 1) and ( 2) in Ω with the piecewise constant conductivity σ , the boundary condition (3) on ∂ Ω, and the initial condition (4) uniquely define the spread of the action potential in Ω. We introduce (u (k,0) , w (k,0) ) as the restriction of the initial data (u 0 , w 0 ) to the layer Ω (k) , and (u (k) , w (k) ) as the restriction of the solution (u, w) to Ω (k) . Each of these restrictions is solution of the monodomain equations set in the subdomain Ω (k) (k = 1, 2):

(x) = σ (1) (x) if x ∈ Ω (1) and σ (x) = σ (2) (x) if x ∈ Ω (
A C∂ t u (k) + f (u (k) , w (k) ) = div x σ (k) ∇ x u (k) + σ (k) 3 ∂ zz u (k) , (7)
∂ t w (k) + g(u (k) , w (k) ) = 0 (8)
for t > 0 and x ∈ Ω (k) with the boundary and transmission conditions

σ (k) ∇u (k) • n = 0 on Γ (k) , k = 1, 2, (9)
σ (1) ∇u (1) • n Σ = σ (2) ∇u (2) • n Σ , u (1) = u (2) on Σ, ( 10
)
where n is the unit normal to ∂ Ω outward of Ω and n Σ is the unit normal to Σ pointing from Ω (2) to Ω (1) , and the initial conditions

u (k) (0, x) = u (k,0) (x), w (k) (0, x) = w (k,0) (x) a.e. x ∈ Ω (k) , k = 1, 2. (11)

Dimensionless monodomain problem

A dimensionless version of the monodomain equations is derived in order to study the balance between the different terms of the equation. We follow the method previously introduced in [START_REF] Colli Franzone | Wavefront propagation in an activation model of the anisotropic cardiac tissue: asymptotic analysis and numerical simulations[END_REF][START_REF] Keener | An eikonal-curvature equation for action potential propagation in myocardium[END_REF][START_REF] Rioux | Numerical Computations of Action Potentials for the Heart-torso Coupling Problem[END_REF]. We first set time and space characteristic lengths t 0 and x 0 , and we use the thickness h > 0 of the two layers to define the dimensionless variables

x = x x 0 , t = t t 0 and z = z h .
For instance, to observe the whole action potential, we will set the length scale to x 0 = 1cm and the time scale to t 0 = 400ms, because they are the typical extent of atrial structures such as pulmonary veins and the duration of an action potential. The two-layers domain

Ω = ω × (-h, h) is mapped to ω × (-1, 1) where ω = 1 x 0 ω ⊂ R 2 .
Using the typical amplitude of an action potential δ u, the resting potential u r , the maximum value of the sodium current f 0 (specifically the maximum value reached by the function f ) and the characteristic conductivity scale σ 0 , we recast the physical quantities, which reads: k) , w (k) ), and ḡ( ū(k) , w(k) ) = t 0 g(u (k) , w (k) ).

ū(k) ( t, x, z) = u (k) (t, x , z) -u r δ u , σ (k) = 1 σ 0 σ (k) , σ (k) 3 = σ (k) 3 σ 0 , f ( ū(k) , w(k) ) = 1 f 0 f (u (
The monodomain equations for the two-layers domain finally read

AC x 2 0 σ 0 t 0 ∂ t ū(k) + f 0 t 0 Cδ u f ( ū(k) , w (k) ) = div x σ (k) ∇ x ū(k) + x 2 0 h 2 σ (k) 3 ∂ zz ū(k) , ∂ t w(k) + ḡ( ū(k) , w (k) ) = 0,
for t > 0 and ( x, z) in the domains Ω(1) := ω × (0, 1) and Ω(2) := ω × (-1, 0). The boundary and transmission conditions ( 9) and ( 10) remain unchanged but are now stated on Γ

(k) = ∂ Ω(k) ∩ ∂ Ω and Σ = ω × {0}.
As a consequence, the dimensionless numbers 1: Electrical conductivities taken from the literature (mScm -1 -monodomain conductivities are the half of the harmonic averages of bidomain ones). Article 1: [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF], article 2: [START_REF] Potse | A comparison of monodomain and bidomain reaction-diffusion models for action potential propagation in the human heart[END_REF], article 3: [START_REF] Clements | Activation dynamics in anisotropic cardiac tissue via decoupling[END_REF] (review), and article 4: [START_REF] Colli Franzone | Wavefront propagation in an activation model of the anisotropic cardiac tissue: asymptotic analysis and numerical simulations[END_REF].

α = AC x 2 0 σ 0 t 0 , β = f 0 t 0 Cδ u and ε = h x 0 Article 1 2 3 4 Intracellular σ 1 3.0 3.0 1.
characterize the solutions. The order of magnitude of each of the physical quantities is set according to usual values for cardiac electrophysiology (see Table 3.2 for some values of the conductivities):

A = 1000cm -1 , C = 1µFcm -2 , δ u = 100mV, f 0 = 100µAcm -2 , σ 0 = 1.5mScm -1 .
The orders of magnitude of the observation scales (t 0 and x 0 ) are then set to x 0 = 1cm, t 0 = 400ms, and consequently, the dimensionless parameters α and β are α = 1.67, and β = 400. Since we are focusing on layers of cardiac tissue, the important quantity is the aspect ratio ε = h/x 0 , assumed to be small. The dimensionless monodomain equations in each layer finally read:

α ∂ t ū(k) + β f ( ū(k) , w (k) ) = div x σ (k) ∇ x ū(k) + σ (k) 3 ε 2 ∂ zz ū(k) , (12)
∂ t w(k) + ḡ( ū(k) , w (k) ) = 0, ( 13
)
and every parameter is fixed and of order 1 except ε and β . According to our analysis, the transverse diffusion ∂ zz ū(k) dominates the reaction terms whenever β ε 2 1. In the next section, we will derive an asymptotic model for ε → 0 while β is kept fixed, consequently in the dominant transverse diffusion regime.

Remark 1. The documented thickness of human atrial tissues is small, but it remains of order h = 0.1cm. Hence, the reaction and diffusion terms actually balance during the depolarization phase (β ε 2 1) and the behavior of the asymptotic model may be different from the observation in atrial tissues. For this reason, we develop a higher order asymptotic model in section 5. During the plateau and repolarization phases, the ionic currents are far less intense, so that the asymptotic regime applies quite well for human atria (β decreases, so that β ε 2 ≤ 1).

Remark 2. Note that the case of tissues with different thickness of the endo and epicardial layers can be dealt with. If we have Ω (1) = ω×]0, λ (1) h[ and k) . We then obtain the equation ( 12), but the transverse diffusion coefficient is now

Ω (2) = ω×] -λ (2) h, 0[, for a given λ (k) > 0, we can define z = z λ (K) h in Ω (
σ (k) 3 (λ (k) ε) 2 and
the current continuity condition through the interface

σ (1) 3 λ (1) ε ∂ z ū(1) = σ (2) 3 λ (2) ε ∂ z ū(2) .
4 Asymptotic one layer model

In this section, we propose a formal asymptotic expansion of the solution of ( 12). The usual surface model of the atria, classical for numerical atrial simulations (see e.g. [START_REF] Haissaguerre | Atrial fibrillatory cycle length: computer simulation and potential clinical importance[END_REF], ROTTER et al., 2007]), is then recovered. Dropping the ¯above the dimensionless quantities, the three-dimensional dimensionless problem reads

α ∂ t u (k) ε + β f u (k) ε , w (k) ε = div x σ (k) ∇ x u (k) ε + σ (k) 3 ε 2 ∂ zz u (k) ε , (14)
∂ t w (k) ε + g u (k) ε , w (k) ε = 0, ( 15
)
for k = 1, 2, t > 0 and (x, z) ∈ Ω (k) . The depence on ε of the solution is emphasized by the subscript. The boundary and transmission conditions read andσ (1)

σ (1) ∇ x u (1) ε • n = 0 in ∂ ω × (0, 1)
3 ∂ z u (1) ε = 0 in ω × {1}, (16) σ (2) ∇ x u (2) ε • n = 0 in ∂ ω × (-1, 0) and σ (2) 3 ∂ z u (2) ε = 0 in ω × {-1}, (17) σ (1) 3 ∂ z u (1) ε = σ (2) 3 ∂ z u (2) ε , u (1) ε = u (2) ε in ω × {0}, (18)
and we assume, for the sake of simplicity, that the initial data are homogeneous in the z direction for each k ∈ {1, 2}, specifically:

u (k) ε (0, x, z) = u 0 (x), w (k) ε (0, x, z) = w 0 (x) in Ω (k) , k = 1, 2, ( 19
)
where the functions u 0 (x) and w 0 (x) are defined in ω.

We then have the following theorem :

Theorem 1 (Coefficients of the asymptotic expansion). If an asymptotic expansion of u

(k)
ε and w

(k)
ε respectively to ε 2 can be exhibited, i.e., if we can write u

(k) ε = u (k) 0 + ε 2 u (k) 1 + ε 4 u (k) 2 + o(ε 4 ), w (k) ε = w (k) 0 + ε 2 w (k) 1 + o(ε 2 ), (20)
then, the asymptotic coefficients (u

(k) 0 , w (k) 0 ) and (u (k) 1 , w (k)
1 ) can be computed explicitly in the following way.

For all t ≥ 0, a.e. x ∈ ω and z ∈] -1, 1[, we have u

(1) andw (1)

0 (t, x, z) = u (2) 0 (t, x, z) := u 0 (t, x)
0 (t, x, z) = w (2) 0 (t, x, z) := w 0 (t, x), where (u 0 , w 0 ) is solution of the surface monodomain problem α(∂ t u 0 + β f (u 0 , w 0 )) = div x (σ m ∇ x u 0 ) , (21)
∂ t w 0 + g(u 0 , w 0 ) = 0, ( 22
)
with the boundary condition σ m ∇ x u 0 • n = 0 on ∂ ω and for t > 0 and initial condition u 0 (0, x) = u 0 (x) and w 0 (0, x) = w 0 (x) in ω where σ m = σ (1) +σ ( 2)

2
is the arithmetic average of the conductivity matrices in both layers.

The function u

(k) 1 is explicitly given by u (k) 1 = b σ (k) 3 z 1 - |z| 2 + c, (23)
with b = div x σ d ∇ x u 0 and c = ū1 -1 6 b σ (2) 3 -σ (1) 3 σ (1)
3 σ

(2) 3

where

σ d = σ (1) -σ (2)
2 and ū1 is solution to the problem

α (∂ t ū1 + β ∇ f (u 0 , w 0 ) • ( ū1 , w1 )) = div x (σ m ∇ x ū1 ) + div x 1 3σ h 3 σ d ∇ x b , (24) ∂ t w1 + ∇g (u 0 , w 0 ) • ( ū1 , w1 ) = 0, ( 25
)
with the boundary condition

σ m ∇ x ū1 • n + 1 2σ h 3 σ d ∇ x b
• n = 0 on ∂ ω and for t > 0, and the initial conditions ū1 (0, x) = 0 and w1 (0, x) = 0 for a.e. x ∈ ω, with

σ h 3 = 2 σ (1) 3 σ (2) 3 σ (1) 3 +σ (2) 3
the harmonic average of the transverse conductivities.

The function w (k)

1 reads:

w (k) 1 (t, x, z) = - t 0 ∂ 1 g (u 0 (s, x), w 0 (s, x)) u (k) 1 (s, x, z) exp - t s ∂ 2 g (u 0 (τ, x), w 0 (τ, x)) dτ ds. (26)
Remark 3. The equations ( 21) and ( 22) are the usual surface monodomaim model for the atria. They are defined on the surface ω and are independent of ε. We note that the surface diffusion tensor is derived by averaging the 3D diffusion tensor in the thickness of the tissue. That surface model has been rigorously derived for instance in [START_REF] Chapelle | A surfacebased electrophysiology model relying on asymptotic analysis and motivated by cardiac atria modeling[END_REF] or in [START_REF] Coudière | An asymptotic two layers monodomain model of cardiac electrophysiology in the atria[END_REF].

Remark 4. We will designate the asymptotic approximation

ũ(k) ε := u 0 + ε 2 u (k) 1 , w(k) ε := w 0 + ε 2 w (k) 1 , (27)
as the enhanced monodomain surface model. This denomination highlights the fact that (

ũ(k) ε , w(k) ε
) is a correction of the usual monodomain surface solution (u 0 , w 0 ) that improves the accuracy of the approximation, as stated in Theorem 2 below.

Proof. We inject the expansions (20) in the system of eqs. ( 14) and ( 15) and we identify the coefficients with the same order, respectively to ε 2 .

First looking at the highest order coefficient, i.e. 1/ε 2 , we obtain the equation

∂ zz u (k) 0 = 0 for k = 1
, 2 together with the boundary and transmission conditions,

∂ z u (1) 0 (t, x, 1) = ∂ z u (2) 0 (t, x, -1) = 0, σ (1) 3 ∂ z u (1) 0 (t, x, 0) = σ (2) 3 ∂ z u (2) 0 (t, x, 0), u (1) 0 (t, x, 0) = u (2) 0 (t, x, 0).
We can easily show that u

(1)

0 = u (2)
0 are constant functions in the z direction. For all z ∈] -1, 1[, we then have u

(1)

0 (t, x, z) = u
(2) 0 (t, x, z) := u 0 (t, x), for t ≥ 0 and a.e.

x ∈ ω.

Using a Taylor expansion of f and g around (u 0 , w 0 ), we get:

f u (k) ε , w (k) ε = f u 0 , w (k) 0 + ε 2 ∇ f u 0 , w (k) 0 • u (k) 1 , w (k) 1 + o ε 2 , g u (k) ε , w (k) ε = g u 0 , w (k) 0 + ε 2 ∇g u 0 , w (k) 0 • u (k) 1 , w (k) 1 + o ε 2 .
The identification of the terms of order ε 0 then leads to the following equation on u

(k) 1 , for k = 1, 2. σ (k) 3 ∂ zz u (k) 1 = α ∂ t u 0 + β f u 0 , w (k) 0 -div x σ (k) ∇ x u 0 , (28)
∂ t w (k) 0 + g u 0 , w (k) 0 = 0, (29)
together with the boundary and interface conditions ( 16) to (18) on u

(k)

1 . The function (t, x, z) → g(u 0 (t, x), .), which is independent of z, completely defines the dynamics of w

(k) 0 . Since, for k = 1, 2, the initial condition w (k) ε (0, x, z) = w 0 (x) is
constant in the z direction through the tissue, the functions w (k) 0 are independent of z for all time t > 0 and have the same value in both layer. We then have w 0 (t, x) := w

(1)

0 (t, x, z) = w
(2) 0 (t, x, z) for all t ≥ 0 and a.e. x ∈ ω, which is solution to equation ( 29). Afterwards, using the boundary conditions ( 16) and ( 17), the integration of equation ( 28) on u (k)

1 along z leads to -σ (1) 3 ∂ z u (1) 1 (., 0) = α (∂ t u 0 + β f (u 0 , w 0 )) -div x σ (1) ∇ x u 0 , (30) σ
(2)

3 ∂ z u (2) 1 (., 0) = α (∂ t u 0 + β f (u 0 , w 0 )) -div x σ (2) ∇ x u 0 . (31)
Adding these two equations with the transmission condition (18) gives the system of equations ( 21)-( 22) on (u 0 , w 0 ), which proves the first part of the theorem.

We then observe that the right-hand side of equations ( 28) is independent of z. Then, the functions z → u (k) 1 (t, x, z) are second degree polynomials. There exists

a (k) , b (k) and c (k) such that u (k) 1 (t, x, z) = a (k) (t, x)z 2 + b (k) (t, x)z + c (k) (t, x).
Introducing Equation ( 21) in ( 28), the right-hand side reads div x (σ m ∇ x u 0 )div x σ (k) ∇ x u 0 which is both equal to ± div x σ d ∇ x u 0 and to the right-hand side of equations ( 30) and ( 31).

Furthermore, the left-hand side of equations ( 28), (30), and (31) are, respectively, 2a (k)

σ (k) 3 , -σ (1) 3 b (1) and σ (2) 3 b (2) . Consequently, we can define b := div x σ d ∇ x u 0 = -2a (1) σ (1) 3 = σ (1) 3 b (1) = 2a (2) σ (2) 3 = σ (2) 3 b (2) .
The potential continuity condition on the interface indicates that we can identify c (1) and c (2) := c. We then have

u (k) 1 = b σ (k) 3 z 1 - |z| 2 + c where b = div x σ d ∇ x u 0 , (32)
and c = c(t, x) is an unknown function. Noting the averaging operator through each layer with •, we have

ū(1) 1 := 1 0 u (1) 1 (•, z)dz = c + 1 3 b σ (1) 3 and ū(2) 1 := 0 -1 u (2) 1 (•, z)dz = c - 1 3 b σ (2) 3
.

As a consequence, we have the relations

ū1 := ū(1) 1 + ū(2) 1 2 = c + 1 6 b σ (2) 3 -σ (1) 3 σ (2) 3 σ (1) 3 , ū(1) 1 - ū(2) 1 2 = 1 6 b σ (2) 3 + σ (1) 3 σ (2) 3 σ (1) 3 := 1 3 b σ h 3 ,
where ū1 denotes the average of u (k)

1 through the whole thickness of the tissue and

σ h 3 = 2 σ (1) 3 σ (2) 3 σ (1)
3 +σ

(2) 3

is the harmonic average of the transverse conductivities. Hence, knowing ū1 uniquely defines c. In the same way, we define w1 as the average of w (k)

1 through the whole thickness of the tissue:

w1 = 1 2 1 0 w (1) 1 (•, z)dz + 0 -1 w (2) 1 (•, z)dz .
Now, in order to close the system on (u

(k) 1 , w (k)
1 ), we want to exhibit a system of equations providing the functions ( ū1 , w1 ).

We then identify the coefficients of order ε 2 in the expansion of (u (k) , w (k) ). We then obtain the equations

σ (k) 3 ∂ zz u (k) 2 = α ∂ t u (k) 1 + β ∇ f (u 0 , w 0 ) • u (k) 1 , w (k) 1 -div x σ (k) ∇ x u (k) 1 , (33) ∂ t w (k) 1 + ∇g (u 0 , w 0 ) • u (k) 1 , w (k) 1 = 0, (34)
with the boundary and transmission conditions ( 16) to (18) on u (k)

2 . We again integrate equation ( 33) for z ∈ (0, 1) and z ∈ (-1, 0), add the resulting equations, use the transmission conditions ( 16) and (17) on u (k) 2 . We then get

1 2 div x σ (1) ∇ x ū(1) 1 + σ (2) ∇ x ū(2) 1 = div x σ m ∇ x ū1 + σ d ∇ x ū(1) 1 - ū(2) 1 2 .
We remark that

ū(1) 1 - ū(2) 1 2 = 1 3 b σ h 3
, and we finally obtain the equations ( 24) and ( 25)

for ( ū1 , w1 ). That gives an explicit computation of u (k)

1 and demonstrates the second part of the theorem.

To end the proof, we observe that the function w

(k)
1 is uniquely determined for k = 1, 2 by equation ( 34) as the solution to a first order linear Cauchy problem of the form

w (t) + a(t)w(t) = -b(t) with a(t) = ∂ 2 g (u 0 (t, x), w 0 (t, x)) and b(t) = ∂ 1 g (u 0 (t, x), w 0 (t, x)) u (k) 1 (t, x, z) and with w(0) = 0 because w (k) 1 (0, x, z) = 0. Its solution is computed explicitly: w(t) = -t 0 b(s) exp -t s a(τ)dτ ds, which is exactly the expression of w (k)
1 given in the Theorem.

5 Asymptotic two layers model 5.1 Averaged equations through the thickness of each layer

In this section, we focus on equations verified by the averages in z ∈ (0, 1) and z ∈ (-1, 0) of the 3D solutions (u

(k) ε , w (k)
ε ) for each of the two layers. Integrating equations ( 33) and ( 34) for z ∈ (0, 1) and z ∈ (-1, 0) we get

α ∂ t ū(k) ε + β f ū(k) ε , w(k) ε = div x σ (k) ∇ x ū(k) ε + (-1) k 3 2 σ h 3 ū(1) ε - ū(2) ε ε 2 + E ε (u), (35)
and

∂ t w(k) ε + g ū(k) ε , w(k) ε = Ē(k) ε (g), (36)
where E ε (u) and Ē(k) ε (g) are remainder terms. In a separate paper [START_REF] Coudière | An asymptotic two layers monodomain model of cardiac electrophysiology in the atria[END_REF] they are shown to be of order ε 3 . The definition of our bilayer model directly results from that observation.

Definition of the two layers model

For given initial conditions û(k) 0 (x) and ŵ(k) 0 (x) on ω, we define the two layers model as the following coupled systems of monodomain equations, for k = 1, 2:

α ∂ t û(k) ε + β f û(k) ε , ŵ(k) ε = div x σ (k) ∇ x û(k) ε + (-1) k 3 2 σ h 3 û(1) ε - û(2) ε ε 2 , (37) ∂ t ŵ(k) ε + g û(k) ε , ŵ(k) ε = 0, ( 38
) with boundary conditions σ (k) ∇ x û(k) • n = 0 on ∂ ω for t > 0 and initial conditions û(k) (0, x) = û(k) 0 (x) and ŵ(k) (0, x) = ŵ(k) 0 (x).
It can be seen as a perturbation of the averaged problems ( 35)-( 36) that neglects the remainder terms E ε (u) and Ē(k) ε (g). We will see in the next section that this perturbation is small enough to correctly approximate the averaged solutions

( ū(k) ε , w(k) ε ) with the bilayer solution ( û(k) ε , ŵ(k) ε )
. Remark 5. We emphasize that, unlike in the enhanced monolayer model, the source function is evaluated in each layer, thus playing a driving role in transmural electrical discrepancies : the bilayer model intrinsically allows inter-layers electrical disconnection and the electrical dynamics can differ from one layer to the other. This is a source of complex propagation patterns.

Error estimate for the surface models

For the completeness of the present article, we introduce two convergence results for both the enhanced monodomain surface model and the bilayer model that we do not demonstrate in that study, but that are fully addressed in a separate paper [START_REF] Coudière | An asymptotic two layers monodomain model of cardiac electrophysiology in the atria[END_REF].

Error estimate on the enhanced monolayer model

We first state that the solutions of the enhanced monodomain surface model ( 27) approximate the 3D solutions with an accuracy of order ε 3 .

We suppose that the initial conditions are constant in the thickness of the tissue. Namely, we assume that u

(k) ε (0, x, z) = u 0 (x) and w (k) ε (0, x, z) = w 0 (x) for (x, z) ∈ Ω (k) , k = 1, 2.
Theorem 2 (Error estimates for the enhanced monolayer model). Supposing that the functions f and g are C 2 (R × R m ) functions and that the asymptotic solutions (u 0 , w 0 ), (u

(k) 1 , w (k) 1 ) (resp. 3D solutions (u (k) ε , w (k) ε ) ε>0
) have a uniform upper bound M respectively to time (resp. to ε and time), we then have, for all 0 < ε ≤ 1, for all T > 0, for 0 ≤ t ≤ T and k = 1, 2, the following estimates

u (k) ε (t) - ũ(k) ε (t) L 2 (Ω (k) ) ≤ ε 3 k 0 , (39)
w (k) ε (t) - w(k) ε (t) (L 2 (Ω (k) )) m ≤ ε 3 k 0 , (40)
∇ x u (k) ε -ũ(k) L 2 (0,t;L 2 (Ω (k) )) ≤ ε 3 k 1 , (41)
∂ z u (k) ε -ũ(k) L 2 (0,t;L 2 (Ω (k) )) ≤ √ 2ε 4 k 1 , (42)
where k 0 and k 1 are independent of ε and only depend on α, σ , Λ 0 , Λ 1 , M, |ω| and T .

Remark 6. We note that Theorem 2 can be seen as the reciprocal of Theorem 1: if we note (u 0 , w 0 ) the solution of the problem (21)-( 22), ( ū1 , w1 ) the solution of (24)-( 25) and (u

(k) 1 , w (k)
1 ) the function defined by ( 23), ( 26), Theorem 2 shows that ( ũ(k) ε , w(k) ε ) are then a third order approximation of the 3D solutions u

(k)
ε and w (k) ε , i.e. that the asymptotic expansion (20) holds, with the norm therein.

Error estimate on the two-layers model

We also state that the two-layers solution ( û(k) ε , ŵ(k) ε ) of ( 37) and ( 38) approximates the averages on each layer ( ū(k) ε , w(k) ε ) of the 3D solutions ( 14) and ( 15) with an order ε 3 accuracy.

We suppose that, for k = 1, 2, we have û(k) 0 (x) = u 0 (x), and ŵ(k) 0 (x) = w 0 (x). Theorem 3 (Error estimates: two-layers model). Assuming that for T > 0 there exists a bound M such that û(k) ε (t, x), ŵ(k) ε (t, x) ≤ M , and under the assumptions of Theorem 2, we have the following estimates: for all 0

< ε ≤ 1, for 0 ≤ t ≤ T , for k = 1, 2, ū(k) ε (t) - û(k) ε (t) L 2 (ω) ≤ ε 3 k 4 , (43)
w(k) ε (t) - ŵ(k) ε (t) [L 2 (ω)] m ≤ ε 3 k 4 , (44)
∇ x ( ū(k) ε - û(k) ε ) L 2 (0,T ;L 2 (ω)) ≤ ε 3 k 5 , (45)
with k 4 and k 5 not depending on ε but only on α, σ , Λ 0 , Λ 1 , M, |ω| and T .

Remark 7. This theorem guarantees that the average by layer in the thickness of the three-dimensional potential converges toward the solution of the two-layers model. Furthermore, the accuracy of the two-layers model is limited by the precision of the approximation of the transverse diffusion which depends in our formulation on the precision of the asymptotic expansion of u (k)

ε . Let us give a version of the problem (37) in physical variables :

A(C∂ t u (k) + f (u (k) , w (k) )) = div x (σ (k) ∇ x u (k) ) + (-1) k γ th (u (1) -u (2) ), (46)
∂ t w (k) + g(u (k) , w (k) ) = 0, (47)
where

γ th = 3 h 2 σ (1) 3 σ (2) 3 σ (1)
3 +σ

(2) 3 is called the coupling coefficient and A, C, σ (k) ,σ

3 are the same parameters as in problem (1), (2).

Remark 8. The model also has a purely electric interpretation, as in [START_REF] Jacquemet | A biophysical model of atrial fibrillation and electrograms : formulation, validation and applications[END_REF]: the two layers are coupled pointwise by a conductance. The coupling coefficient from our model, γ th , is expressed in mS cm -3 , which is consistent. Consequently, the term γ th (u (1)u (2) ) is a density of current expressed in µA cm -3 . Remark 9. We recall that this approximation is expected to be correct if the thickness of the atrial tissue is of order 0.01cm, which is not a realistic range for human hearts. Indeed, surface atrial models usually takes into account only the first order two-dimensional monodomain approximation (u 0 , w 0 ), and cannot capture complex propagations patterns like observed in figure (3) below.

Remark 10. Following Remark 2, we can perform the same computation to derive the asymptotic coupling coefficient when Ω (1) and Ω (2) have respectively a thickness λ (1) h and λ (2) h, with λ (h) > 0 and of order 1. Equations ( 46) hold with

γ (k) th = 3 2 σ h λ λ (k) h 2 , where σ h λ = 2 λ (1) σ (1) 3 + λ (2) σ (2) 3 -1
. We note that σ h λ is the harmonic mean of the transverse diffusion coefficient weighted by the thickness coefficient and that the coupling coefficient γ

(k)
th differs on each layer due to the factor λ (k) .

Numerical illustrations

In order to complete the previous theoretical foundations, we focus on the potentatiality of the two-layers models in reproducing some actual complex 3D phenomena in a physiological context. This is a threefold objective: analysis of the asymptotic regime of the surface models with respect to the physiological thickness of atrial tissues; introduction of an optimization heuristic so as to correct the behaviour of the model for physiological thicknesses; study of the robustness of the model for physiologically relevant 3D geometries.

Numerical method

Test cases: general protocol

The general protocol is as follows. First we compute the solution (u

(k) 3D , w (k)
3D ) of the 3D problem (4) in a domain Ω = Ω (1) ∪ Ω (2) with a diffusion tensor σ discontinuous across the interface Σ. Namely, the diffusion tensor σ is defined as in equations ( 6) and ( 5), and we choose the fibres to be perpendicular one with another: θ (2)θ (1) = π 2 . This solution serves as a reference. Then, we compute the solutions to the three following surface models: 1×2D: the solution (u 0 , w 0 ) of the classical monolayer surface model ( 21) -( 22)

with diffusion tensor σ m = σ (1) +σ (2)

2

. In this case, the solution (u 0 , w 0 ) is isotropic;

2×2D-th: the solution ( û(k) γ th , ŵ(k) γ th )
of the two-layers model ( 46) -( 47) with the theoretical coupling coefficient γ th given in the text;

2×2D-opt: the solution ( û(k) γ opt , ŵ(k) γ opt )
of the two-layers model ( 46) -( 47), but with a coupling coefficient γ opt defined in a physiologically optimal manner, see section 7.1.5. Indeed, we will observe that the 1×2D and 2×2D-th models are no longer accurate for physiological thicknesses of atrial tissues. Hence we introduce the third, optimized model, so as to improve the accuracy of the 2×2D-th for thicknesses within the physiological range (that is, outside of the asymptotic regime).

The reference geometry

We use a reference perpendicular geometry, as defined in the previous theoretical study, in order to investigate the convergence and accuracy of the two-layers model. More general geometries will be described in section 7.1.6. Slab of tissue. The 3D slab is the domain Ω = ω × (-h, h) with ω = (0, L x ) × (0, L x ), of length L x = 1cm, and half-thickness h varying between 0.005 cm and 0.2 cm. The total thickness, 2h, compares to physiological scales varying from the length of one cardiomyocyte (100µm) to the average thickness of human atria (about 0.4 cm). The fibres directions are set to ν

(1) 1 = (1, 0, 0) and ν

(2) 1 = (0, 1, 0) (see equation ( 5)).

We consider a cartesian grid on Ω (resp. on ω) with discretization steps ∆x × ∆y×∆z (resp. ∆x×∆y). The 3D and 2D grids are decomposed in simplicial meshes (i.e. with respectively tetrahedral and triangular elements). Each voxel of the 3D grid is decomposed in 24 tetrahedra whereas the voxels of the 2D grid are decomposed in 4 triangles, as presented in figure 7.1.2. We take ∆x = ∆y = 5.10 -3 cm, which is small to guarantee the numerical convergence of the discrete solution. To check the convergence of the 2×2D model when h vanishes, we vary h and ∆z following two constraints. We impose that each mesh contains at least 10 voxels in the thickness, in order to capture any transversal activity, even for small thicknesses. For large thicknesses, we keep the number of degrees of freedom in a tractable level while preserving the mesh element shapes, by taking ∆z ≤ 4∆x. The final value of ∆z for each value of h are presented in Table 3. We denote by DOF3D (resp. DOF2D) the number of degrees of freedom in Ω (resp. ω).

Electrophysiological background

We chose the Beeler-Reuter ventricular model [START_REF] Beeler | Reconstruction of the action potential of ventricular myocardial fibres[END_REF], which is computationally simple but reproduces correctly the depolarization dynamics of cardiomyocytes. All the other parameter values are given in Table 7.1.3.

A stimulation was applied by fixing the potential to V m,init = 20mV for 0 ≤ t ≤ 2ms, in the square S = 1 2 (L x , L x ) + (-2∆x, +2∆x) × (-2∆x, +2∆x) for the 2D test cases, and 1 2 (L x , L x , 0) + (-2∆x, +2∆x) × (-2∆x, +2∆x) × (-h, +h) in the 3D test cases.

Resolution

All the models are discretized by the standard P1-Lagrange finite element method with mass lumping and the Rush-Larsen time-stepping scheme [START_REF] Rush | A practical algorithm for solving dynamic membrane equations[END_REF]. We used a fixed time-step ∆t = 0.05ms. The diffusion terms are discretized implicitely, and the resulting linear system is solved with a Jacobi preconditionned conjugate gradient algorithm with a fixed tolerance equal to 10 -10 .

Optimizing the coupling coefficient

According to Theorem 3, the difference between the two-layers and the 3D models decreases in O(ε 3 ) when the coupling coefficient is γ th = 3 2 σ h 3 ε 2 . Anyway, for most physiological thicknesses, we will observe that this difference remains qualitatively too high. Hence, for such a thickness (fixed), we modify the coupling coefficient in order to minimize this difference, resulting in the 2×2D-opt model.

To this aim, we don't use the functional norms from Theorem 3, but the more relevant activation maps. Activation maps are the spatial disctributions of the timeof-arrival of the activation wavefront, which are a major indicator in clinical and biological practice. For a given threshold θ ∈ R, and a solution u

∈ C 0 ([0, T ], L 2 (ω)), the activation map T θ (u) ∈ L 2 (ω) is defined by T θ (u)(x) = inf {t ∈ (0, T ), s.t. u(t, x) > θ }.
Note that the activation maps are defined on the surface ω. In practice, we define one activation map for each layer, defined for θ = -40mV.

For given 3D solutions u (k)

ε , k = 1, 2, we define the cost function to minimize by σ iso cm -1 µF cm -2 mS cm -1 mS cm -1 500 1 1.5 0.2 0.2 0.85 46) -( 47) computed with the coupling coefficient γ, and ū(k) 3d is the 3D solution averaged in the thickness of the layer k.

J(γ) := T θ ū(1) 3d -T θ û(1) γ L 2 (ω) T θ ū(1) 3d L 2 (ω) + T θ ū(2) 3d -T θ û(2) γ L 2 (ω) T θ ū(2) 3d L 2 (ω) , thickness (
We finally look for a value of γ in the neighbourhood of the theoretical one, γ th , that achieves a local minimum of J. Therefore, we use the Brent's method implemented in the scipy library [START_REF] Oliphant | Python for scientific computing[END_REF], that is interfaced with our solver.

Remark 11 (Gold standard). We note that the design of the 3D gold standard is a key issue of the optimization process. In order to construct a physiological 2×2D model of the atria, the 3D gold standard should: 1) be simple to compute, and 2) should contain all the relevant transmural heterogeneities that we want to capture. The 3D slab previously defined is a good candidate for such a gold standard. However, there is a priori no reason why the optimized coupling coefficient obtained with this simplified 3D reference should give satisfactory results on the more complex geometries that are observed in the atria. We then check the robustness of that optimization strategy on physiologically relevant additional geometries: a curved geometry, and a geometry with a non homogeneous thickness. We also verify if the optimization method is improved when the 3D gold standard is computed directly on the complex 3D geometry, instead of the simplified slab of tissue.

Additional geometries

These geometries will be used to challenge the two-layers model in more relevant situations.

Curved geometry. It is the domain Ω = Ω (1) ∪ Ω (2) , where Ω (k) are defined in cylindrical coordinates by,

Ω (1) = {(r, ψ, z), R -h < r < R, 0 ≤ ψ ≤ 2π, 0 < z < 1}, and Ω (2) = {(r, ψ, z), R < r < R + h, 0 ≤ ψ ≤ 2π, 0 < z < 1}. Hence, the interface is now the surface ω = {(R, ψ, z), 0 ≤ ψ ≤ 2π, 0 < z < 1}.
The fibre directions are choosen to be longitudinal in the layer 1 (ν (1) 1 = ν z ) and circumferencial in the layer 2 (ν

(2) 1 = ν ψ ).
We take R = 1 2π + hcm and h = 0.05cm, so that the domain Ω has the typical dimensions of a human pulmonary vein with basal ring of thickness 0.1cm and minimal perimeter of 1cm, and realistic fibre orientations.

Variable thickness. This domain is a prism with a trapezoid basis It can be defined by Ω = (0, 1) × T in cartesian coordinates, where T = {(y, z), 0 ≤ z ≤ (0.05 + 0.10y), 0 ≤ y ≤ 1}. In this case, the interface is the surface ω = (0, 1) × {(y, z = 0.5 (0.05 + 0.10y)), 0 ≤ y ≤ 1}, and the fibre direction is in the x-direction in Ω (1) (ν (1) 1 = (1, 0, 0)) and along the thickness variation gradient in Ω (1) (ν

(1) 1 = (0, 1, 0.10z)).
Discretization. Both geometries Ω, and both interfaces ω are discretized by a cartesian grid in the parametric space, of size 200 × 200 × 10, which is in turn decomposed with a simplicial mesh. The spatial resolution of these meshes is then 5. 10 -3 cm.

results

We first focus on the qualitative behaviour of the different models. We assess the propagation pattern in the 3D model for different regimes of tissue thickness and observe the corresponding results given by the surface models. We then compare the accuracy of the surface models, and finally study the optimal coupling coefficient in significant physiological situations.

Qualitative behaviour

The 3D model. We first aim at illustrating the 3D propagation pattern in three regimes of tissue thickness: 1) 2h = 0.01cm, which represents a non-physiological thickness of tissue -two myocytes in the thickness of the tissue -but lays in the convergence regime of both surface models -cf section 3.2 -, 2) 2h = 0.1cm, which is the characteristic thickness of thin cardiac tissues like pulmonary veins [START_REF] Ho | Architecture of the pulmonary veins: relevance to radiofrequency ablation[END_REF] and 3) 2h = 0.4cm which is in the range of the documented thickness of atrial tissues [START_REF] Angel Cabrera | The architecture of the left lateral atrial wall: a particular anatomic region with implications for ablation of atrial fibrillation[END_REF].

As expected from Theorem 2, we observe in Figure 2 a uniform potential distribution through the thickness of the geometry for the thin tissue 2h = 0.01cm. On the contrary, for thicker tissues, the propagation follows a very different behaviour. For 2h = 0.4cm, the tissue seems to be thick enough, and the current strong enough, for the wave propagation to be mainly directed towards the fibre direction in each layer. A boundary layer located at the interface of the fibre layers connects that two unidirectional anisotropic waves without significantly perturbing their propagation near the upper and lower boundaries of the slab. For 2h = 0.1cm, the propagation pattern is intermediate: the thickness of the boundary layer connecting the waves of both tissue layers is in the range of the thickness of the slab. It then strongly Following the theoretical results for ε going to zero, the difference between the averaged potentials in each layer should be of order ε 2 at any position x ∈ ω. From the 3D simulations it appears that this is not the case eveywhere. In fact due to the very strong depolarisation current included in the electrophysiological model, very sharp activation fronts appear and are leaning to propagate in a different manner in each layer although the coupling between layers decreases their difference of potential. Around the activation front we observe different but slow varying values of the potential in the thickness of each layer. But near the surface z = 0 these values are joined by a fast varying boundary layer. This discrepancy with the theoretical asymptotic analysis is due to the strong activation current which implies (cf 1) that the assumptions for this analysis are not verified close to the front.

The 2×2D-th and 1×2D models. We now investigate the qualitative behaviour of the 1×2D and the 2×2D-th models for the same three thicknesses. Figure 3 displays a snapshot of each model at the same time t = 11ms for these three thicknesses, with layer 1 on the top row and layer 2 on the bottom one. The 3D solution is averaged in the thickness of each layer before being displayed, while the 1×2D model has only one layer. The isoline V m = -40mV of each model at t = 11ms is also plotted.

For the thin tissue, 2h = 0.01cm, the three models have a very similar quasiisotropic behavior, as expected after Theorems 2 and 3.

For the moderately thick tissue 2h = 0.1cm, the 2×2D-th model reproduces qualitatively well the diamond-shape of the depolarisation front from the 3D model. Though, there remain some important discrepancies, notably in the direction perpendicular to the fibres. The isotropic propagation of the 1×2D model is qualitatively far from the 3D solution.

In the thick tissue, 2h = 0.4cm, the same comments hold: discrepancy in the direction transverse to the fibre for the 2×2D-th model, and irrelevance of the 1×2D model.

As a conclusion, for physiological thicknesses, the theoretical coupling coefficient provides insufficient coupling between the layers: the 2×2D-th model does not reproduce the transmural perturbations of the propagation due to the boundary layers, and therefore the fronts in each layer tend to propagate independently one from another.

Remark 12 (Intermediate potential zones). The averages used to compare the 3D model with the surface models keep track of the boundary layer at the interface. for instance, in Figure 3 The 2×2D-opt does not provide noticeable additional information in the 2h = 0.01cm case, as compared to the 2×2D-th model. For thicker tissues, the 2×2D-th model is improved, mainly in the transverse direction: the comparison between Figures 3 and4 shows that the 2×2D-opt localizes more accurately the front than the 2×2D-th model, even for large thicknesses. It also better reproduces the electrical state of the boundary layer of the 3D model.

Convergence of the surface models

10 -2 10 -1 2h (cm) 10 -2 10 -1
Relative error on activation maps 2x2D-opt 2x2D-th 1x2D

Figure 5: Convergence of the relative L 2 error on activation maps. The relative error measured in L 2 norm on the activation map for the 2×2D-opt (blue), 2×2Dth (green) and 1×2D (red) models vs the thickness of the tissue, in log-log scale.

Figure 5 displays the L 2 (ω) relative error on the activation maps, namely the cost function J for the 2×2D-opt, 2×2D-th and 1×2D models.It shows that the three surface models actually converge towards the 3D solution when h decreases. For the smallest thickness, 2h = 0.01cm, a slight discrepancy persists between the 1×2D and 2×2D-th models: the former has a relative error around 2%, whereas the latter has an error less than 0.8%. Although the difference was not noticeable on Figures 3 and4, the 2×2D-th model is more accurate even for very thin tissues, in accordance with Theorems 2 and 3.

For 2h ≤ 0.02cm, the relative error of the 2×2D-th increases. This is most probably a numerical issue, related to the anisotropy introduced in the 3D meshes for thin layers. The optimization then unduly compensates the numerical errors. But there is no contradiction with Theorem 3, which is stated with L 2 (ω) norm on potential maps, and not on activation maps.

We recall that the 2×2D-th model has a better theoretical order of convergence than the 1×2D model (Theorems 2 and 3). This is confirmed in Figure 5. The slow convergence of the 1×2D model can be furthermore explained: whereas a quasiisotropic propagation is observed in the 3D model for 2h ≤ 0.06, which should promote the convergence of the 1×2D model, the different transient regimes, that immediately follow the stimulation, yields a delay of a few milliseconds between the activation fronts (time-asymptotics).

For 2h ≥ 0.1cm, that is within the physiological range, the relative errors of the 1×2D and the 2×2D-th models are comparable and both superior to 10%. The large error of the 2×2D-th model mainly comes from the bad approximation of the transverse electrical current. The 2×2D-opt model better accounts for the transverse electrical activity, and its error remains under 4.5%. This illustrates the relevance of the optimization procedure proposed in section 7.1.5.

Table 3 finally shows that the computational cost is reduced by a factor in the range 25 to 70 for the 2×2D models, and 70 to 170 for the 1×2D model, with respect to the 3D model.

Study of the locally optimal coupling coefficient

In order to investigate the capabilities of the 2×2D-opt model, we first want to investigate the discrepancies of coupling coefficient between the 2×2D-th and 2×2D-opt models for a given geometry. We next want to check how accurate is the optimal coupling coefficient computed on a slab of tissue in different physiologically relevant geometrical situations. γ th (mS cm -3 ): theoretical coupling coefficient computed with the formula defined in Equation ( 46). γ opt (mS cm -3 ) : optimized coupling coefficient.

Optimal coupling coefficient in a slab of tissue. As observed on Figure 6, for 2h ≤ 0.07cm, the optimal and theoretical coupling coefficients are almost the same. Coupling coefficient (mS.cm

-2 ) γ th γ opt prisme cylinder

Figure 6: Optimized γ opt (green) and theoretical γ (blue) coupling coefficients versus the thickness 2h. The blue and pink bullets represent the optimal coefficient γ opt for the trapezoïdal and cylindrical geometries for a thickness 2h = 0.1cm.

The slight divergence for 2h < 0.03cm may be due to discretization issues as discussed below. For physiological tissue thickness, 2h > 0.1cm, the theoretical coefficient γ th obviously underestimates the coupling between layers, as explained in §7.2.1. Specifically, the theoretical two-layers model underestimates the transverse component of the current in the boundary layers near the interface. This current is exactly the coupling term in the model. Hence the transverse, coupling current, increases linearly with the coupling coefficient. This explains the success of the optimization procedure proposed above.

Robustness in a cylindrical domain. On the cylindrical domain, of thickness 2h = 0.1cm, we can either use the optimal coupling coefficient γ opt found on the parallelepipedic slab of tissue, or compute the optimal coupling coefficient γ opt-cyl suited to this geometry. But at the expense of a new optimization procedure. As observed on Figure 6 (pink bullet), the two coefficient are finally close one to another: γ opt-cyl = 198.9mS cm -3 while γ opt = 185.4mS cm -3 . This justifies that we can use the optimized coefficient from the simpler parallelepipedic slab on a more general geometry. with, the optimized coefficient γ opt , the optimized coefficient γ opt-cyl , and to the 3D model. The corresponding isolines are shown on Figure 8. Both figures prove a very good accordance between the three models. In details, the 2×2D model with γ opt-cyl concentrates the error in the axial direction of the cylinder, while the 2×2D model with γ opt model has its error in the circumferential direction. The relative error on activation maps respectively to the 3D model are precisely: 2.4% for the γ opt-cyl model, and 2.5% for the simpler γ opt model. As a comparison, the relative error found for the γ opt model on the 0.1 cm thick parallelepipedic slab of tissue is equal to 3.3%.

Robustness for variable thickness. Again, we can compute the optimal coupling coefficient γ opt-trap suited to the variable thickness geometry (blue bullet on Figure 6). We find that γ opt-trap = 187.6mS cm -3 , while γ opt = 185.4mS cm -3 has not changed. Note that 2h = 0.1cm is only the mean thickness of the trapezoid geometry. Hence a third model is possible, by assigning a new coupling coefficient γ opt-nh (x, y) depending on the thickness of the tissue at location (x, y). This coefficient can be read on the green line on Figure 6, and is actually computed by linear interpolation from the results of Table 5. Figures 9 and 10 display the snapshots and isolines at t = 11ms for this test case. The homogeneous models, coefficients γ opt = 185.4mS cm -3 (red isoline) and γ opt-trap = 187.6mS cm -3 (green isoline), are indistinguishable. The associated relative error on the activation map is exactly 5.7% for both models. The model with variable coefficient γ opt-nh (x, y) slightly improves the results: it notably reproduces the dissymetry of the propagation, and always improves the localizaton of the front except in the lower part of layer 2. The relative error on the activation maps is equal to 5.1%. Again, the optimized coefficient from the simpler perpendicular slab can be used with reasonable accuracy.

Discussion

The idea of a multi-layer model applied to atrial electrophysiology modelling is not new in itself: a similar two-layers model was proposed in [START_REF] Jacquemet | A biophysical model of atrial fibrillation and electrograms : formulation, validation and applications[END_REF][START_REF] Gharaviri | A computer model of endo-epicardial electrical dissociation and transmural conduction during atrial fibrillation[END_REF]. But its rigorous mathematical derivation, including dimensionless formulation, asymptotic derivation and convergence properties, together with the numerical exploration of its properties are original. The proofs of the convergence theorem can be found in a more theoretical paper [START_REF] Coudière | An asymptotic two layers monodomain model of cardiac electrophysiology in the atria[END_REF] or in [?].

Advantages of the 2×2D models

Mathematical foundations. Both convergence rate and asymptotic behaviour of the 2×2D model are specified. That ensures the consistency of the model and gives us a solid foundation for linking the physical properties of the atrial tissue to the parameter values of the 2×2D models. It also better characterizes the modelling assumption, specially the transmural heterogeneities, then giving a better control on the simplification introduced in the model.

2×2D vs 1×2D models. The 2×2D model supplies several improvements to the usual 1×2D surface models. It gives a better approximation order than the monolayer model, both theoretically -cf Theorems 2 and 3-and numerically -cf Figure 5. It then improves the approximation of the 3D model for small thickness. But overall, the optimization process provides acceptable approximations for physiological thickness of tissue, outside the range of convergence of the surface models.

The key improvement of the bilayer model over the usual surface model lies in the distinction of individual reaction and diffusion terms in the layers. This is important to trigger transmural gradients and dissociation of the electrical activity, together with complex anisotropic propagation patterns as observed in [START_REF] Vetter | Epicardial fiber organization in swine right ventricle and its impact on propagation[END_REF]. The usual surface model, even enhanced by addition of the second order term ε 2 u (k) 1 would not easily account for these phenomena.

2×2D vs 3D models. The main advantage of the 2×2D models is the weak computational load, comparatively to the 3D model. With the numerical protocol applied in this paper, the speed-up is 50 on average. This amount may be overestimated comparatively to other discretizations, due to the choice of simplicial elements in the finite element method. That leads to matrices with larger bandwith than given by other elements, such as quadrangles. Nevertheless, modelling the three dimensional features with a 3D model requires to include a sufficient number of three dimensional elements in the thickness of the tissue, weighting the computational load.

Existing 3D physiological models of the atria are implemented with 3D hexahedral elements [START_REF] Tobón | A three-dimensional human atrial model with fiber orientation. electrograms and arrhythmic activation patterns relationship[END_REF]. Despite the presence of several layers of 3D elements, this kind of model may fail to capture transmural electrical behaviours, due to the small number of transmural freedom degrees. Two alternative approaches can be developed. 1) A sufficient number of freedom degrees can be added through the thickness of the tissue, deteriorating computing speed. 2) The transmural conduction coefficient can be tuned ad hoc. A too small number of transmural elements leads to homogenize the endo-epi distribution of potential. In order to model transmural effects, this homogenisation process can be counterbalanced by an heuristic method similar to the process presented in §7.1.5. The transmural conduction coefficient of the 3D elements could be optimized to fit a 3D propagation computed with a finer transmural space step.

Physiological relevance

Understanding the underlying physics. The dimensional analysis from section 3.1 gives a good understanding of how the trade-off between reaction and diffusion could trigger complex three-dimensional propagation patterns. If the transverse diffusion is predominant, the overall propagation front is isotropic. On the contrary, the propagation front will be diamond shaped when the transverse diffusion and the transmembrane current have a balanced magnitude. As the source function intensity varies through the action potential, the relative strength of the source and diffusion terms also varies. That lead to different front shapes at depolarisation and repolarisation time. The 2×2D models correctly model this phenomenon.

Coupling coefficient. The coupling coefficient gathers several physiological informations such as the transmural conductivity or the thickness of the tissue. The definition of its value is then a key point. This paper gives a theoretical value for vanishing thickness, and a heuristic method, based on an optimization procedure, that enhances the value for physiological thickness -see §7.2.3.

We point out that the choice of the threshold value θ -cf §7.1.5 -that determines the location of both fronts in the optimization process is important. Tuning θ allows us to approximate the pulling effect of the boundary layer which, by connecting the propagation front of both layers, increases the velocity of the propagation orthogonally to the fibre direction in each layer-cf Figure 11 for illustrations.

Physiological models. Endo-epi electrophysiological heterogeneities can easily be included in a 2×2D model.

First, transmural variations of the structure of the tissue such as fibre direction, conduction properties, pathological heterogeneities -fibrosis, ...-are modelled by the layer-dependent conductivity tensor σ (k) . But functional heterogeneities can also be considered through layer-dependent source functions f (k) and g (k) . Likewise, it is possible to deal with atrial tissue where the distinct layers have different thicknesses.

Those transmural heterogeneities can be easily dealt with in the theoretical study of the 2×2D model. Remark 2 and 10 show that a simple rescaling can be included in the expression of the theoretical coupling coefficient to take into account endo-epi thickness difference in the asymptotic regime. But they can also be included in the heuristic optimization of the coupling coefficients. The 3D gold standard must contain the structural or functional properties that we want to model. The 2×2D-opt model is then designed to best fit the 3D propagation. The 2×2D model can also deal with full endo-epi electrical disconnections. That situation is observed e.g. in the pectinate muscle zone of the right atrium: those muscles are attached to the epicardial layer in their extremities but can be free of connections in their middle. The coupling coefficient has just to be set to zero in that region.

The experiments presented in §7.2.3 show that the optimal coupling coefficient γ opt computed with a slab of tissue as gold standard leads to a 2×2D-opt model that can correctly approximate the depolarization spread in a curved or in a thickness-varying geometry. Computing the optimal coupling coefficient with a curved or a thickness-varying 3D gold standard gives very similar results. Moreover, a non homogeneous distribution of the coupling coefficient in the 2×2D-opt model enhances the approximation of the propagation in a trapezoidal geometry. Then, choosing a slab of tissue as the optimization gold standard presents several advantages. First, this method is obviously more tractable than constructing an accurate 3D model of the atria as gold standard. Secondly, this 2×2D-opt model is robust respectively to the main geometrical deformations observed in the atria: curvature and thickness variations.

Limitations

The accuracy of the two-layers model may be enhanced. We could e.g. perform a better approximation of the flow of current through the interface by adding other layers or including the next terms of the asymptotic expansion of the one-layer model.

We chose to work with the monodomain model all along this paper. This choice has been guided by the sake of clarity of the presentation and notations. The twolayers bidomain model can be derived by using the same process, once having defined the flux transmission condition through the interface of both layers, for the intra and extra cellular media.

The transmural heterogeneities that are well rendered by the two-layers model are mainly located in small parts of the atria (pulmonary veins, left atrium posterior and anterior walls, pectinate muscles, and Crista terminalis). Hence, coupling the two-layers model in that zones to the usual surface model in the remaining tissue might improve the computational efficiency without degrading the accuracy of the overall method.

Complicated fibre structures observed in the pulmonary veins, such as blending or crossing layers or tongues of tissue may fail to be described by a two-layers representation. A specific description should be designed to deal with those structures. Furthermore, several transatrial structures such as the Bachmann's bundle may lead to a non manifold geometry in a surface representation. An asymptotic analysis could be performed on a 3D description of the insertion of that bundle in the fibre distribution of the atrial chamber to derive the asymptotic transmission conditions that should be imposed at the interface between the 2D representations of the bundle and the chambers.

Conclusion

We studied the mathematical derivation of a two-layers model of atrial electrophysiology, and we explored numerically its convergence for vanishing thickness. The two-layer model is improved by an optimal tuning of its coupling coefficient. The robustness of this optimal two-layers model to geometrical variations -curvature, non homogeneous thickness-is finally assessed. We then obtain solid mathematical basis to promote the use of such a model in physiological and pathological situations.

This 2×2D model is interesting to construct anatomically and structurally accurate models of the human atria. Such a model should be based on clinical images that give the two-dimensional surfaces, that are complemented by histological descriptions. An optimal coupling coefficient distribution should be constructed, relying on an a priori knowledge of the structure and the function of the atrial tissue.

The 2×2D model gives us a set of tools able to model pathological patterns such as epi-endo electrical dissociation, propagation blocks in one layer or the alteration of the propagation velocity given by endo-epi heterogeneities, which are expected to play a major role in the initiation and the perpetuation of arrhythmias. A first model for that purpose was presented in [START_REF] Labarthe | A bilayer model of human atria: mathematical background, construction and assessment[END_REF]. We believe that it can be an appropriate tool to study the structure-to-function relations in the atria, from a medical point of view.

Figure 1 :

 1 Figure 1: Left: 2D mesh. Middle : 3D mesh related to the 2D one. Right : slice of the discretization of a voxel. A voxel is composed of 24 tetrahedra.

Figure 2 :

 2 Figure 2: Snapshot for t = 11ms of the potential map of the 3D model for different thicknesses 2h. In each view, we clipped the subdomain [0; 0.5] 2 × [0; 2h].

  for 2h = 0.4cm, the potential rises up to an intermediate 2h

Figure 3 :Figure 4 :

 34 Figure 3: Snapshot at t = 11ms for different thicknesses of the transpotential maps for the 2×2D-th, 3D and 1×2D models. Columns (a)-(b): 2×2D-th and 3D models. Top: layer 1. Bottom: layer 2. Column (c): 1×2D model. Column (d): isoline V m = -40mV for the 2×2D-th (blue), 3D (black) and 1×2D model (purple). We recall that the 3D model is averaged through the tickness of each layer, i.e. for z ∈]0, h[ or z ∈]h, 0[.

Figure 7 :Figure 8 :

 78 Figure 7: Snapshot of a depolarization sequence in the 3D and two-layers cylindrical models. Potential maps taken 11ms after stimulation. (a)-(b) : two-layers model with an optimized coupling respectively computed in the slab (γ opt = 185.4mS cm -3 ) and in the cylindrical geometry(γ opt-cyl = 198.9mS cm -3 ). (c) : 3D model.

Figure 9 :Figure 10 :

 910 Figure 9: Snapshot of a depolarization sequence in the 3D and two-layers models. Potential maps taken 11ms after stimulation. (a)-(b) : two-layers model with an optimized coupling respectively in the slab (γ opt = 185.4mS cm -3 ) and in the prism geometry(γ opt-trap = 187.6mS cm -3 ). (c): 2×2D model with a non homogeneous coupling coefficient: in each point, the coupling coefficient is an approximation of the optimal coupling coefficient for the thickness of the tissue in that point. (d): 3D model.

Figure 11 :

 11 Figure11: Significance of the choice of θ . Snapshot of the 3D solution averaged in the thickness of the layer 1 at t = 11ms for 2h = 0.4cm and isolines V m = 0mV (white), -40mV (green) and -70mV (black). As emphasized in Remark 12, the averaging of the 3D potential maps through the thickness of the tissue leads to intermediate potential zones, which reflects the boundary layer potential variation. If θ is too high, e.g. θ = 0mV, the isoline V m = θ does not cross the intermediate potential zone. That would mean that the propagation front that we want to fit does not involve the boundary layer potential. If θ is too low, e.g. θ = -70mV, all the intermediate potential zone is considered as a depolarized tissue, which is not physiological. For θ = -40mV, the isolines include a part of that boundary layer trace and the 2×2D-opt model can take it into account.

Table 3 :

 3

		2h) characteristics	3D	2D bilayer 2D monolayer
		∆z	0.02	-	-
	0.4	# DOF	3 296 421	160 802	80 401
		CPU time ratio	1	0.0171	0.0069
		∆z	0.015	-	-
	0.3	# DOF	3 296 421	160 802	80 401
		CPU time ratio	1	0.0162	0.0066
		∆z	0.01	-	-
	0.2	# DOF	3 296 421	160 802	80 401
		CPU time ratio	1	0.0146	0.0059
		∆z	0.01	-	-
	0.1	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0310	0.0127
		∆z	0.009	-	-
	0.09	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0265	0.0118
		∆z	0.008	-	-
	0.08	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0271	0.0119
		∆z	0.007	-	-
	0.07	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0321	0.0142
		∆z	0.006	-	-
	0.06	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0324	0.0142
		∆z	0.005	-	-
	0.05	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0347	0.0145
		∆z	0.004	-	-
	0.04	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0365	0.0152
		∆z	0.003	-	-
	0.03	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0373	0.0155
		∆z	0.002	-	-
	0.02	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0313	0.0130
		∆z	0.001	-	-
	0.01	# DOF	1 688 411	160 802	80 401
		CPU time ratio	1	0.0216	0.0094

Presentation for each thickness 2h of the characteristics of the discrete domains : length, width, thickness, space steps in each direction (All that dimensions are expressed in cm), number of degrees of freedom and the ratio of the computational time over the three-dimensional one : CPU time / CPU time 3D.

Table 4

 4

		: Physiological parameters used for the different simulations
	where	û(k)

γ , k = 1, 2 are the solutions of the bilayer problem (

Table 5 :

 5 Optimized coupling coefficient for a tissue panel. 2h (cm) : total thickness.

	2h (cm)	0.4	0.3	0.2	0.1	0.09	0.08	0.07
	γ th (mS cm -3 )	7.5	13.3	30	120		148.1	187.5	244.9
	γ opt (mS cm -3 )	68.2	84.4	113.5	185.4	202.3	225.4	261.6
			0.06	0.05	0.04		0.03	0.02	0.01
			333.3	480	750	1333.3 3000 12 000
			326.0	440.5	818.9 1329.9 2479.4 7487.3

Acknowledgements

We used the experimental testbed PLAFRIM , being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr) to perform the computations presented in the article.

This work was partially supported by an ANR grant part of "Investissements d'Avenir" program with reference ANR-10-IAHU-04 and ANR-13-MONU-0004-01, and by the Region Aquitaine.