An asymptotic two-layer monodomain model of cardiac electrophysiology in the atria: derivation and convergence *

. We investigate a dimensional reduction problem of a reaction-diﬀusion system related to cardiac electrophysiology modeling in the atria. The atrial tissues are very thin. The physical problem is then routinely stated on a two-dimensional manifold. However, some electrophysiological heterogeneities are located through the thickness of the tissue. Despite their biomedical signiﬁcance, the usual dimensional reduction techniques tend to average and erase their inﬂuence on the two-dimensional propagation. We introduce a two-dimensional model with two coupled superimposed layers that allows us to take into account three-dimensional phenomena, but retains a reasonable computational cost. We present its mathematical derivation, show its convergence toward the three-dimensional model, and check numerically its convergence speed.

Atrial tissue is very thin.In order to take advantage of this anatomical specificity, atrial electrophysiological models are often formulated by setting the monodomain problem on bidimensional manifolds [5,7].Those surface models can be derived by asymptotic analysis which results in averaging structural and functional informations through the thickness of the tissue; see Remark 2 below or [2,4].However, the transmural distribution of those anatomical and electrical characteristics can be very heterogeneous, despite the thinness of the tissue.For example, abrupt fiber direction discontinuities through the wall are documented [8,16].That can trigger complex propagation patterns [20] which are suspected to be a substrate for arrhythmia [12].In order to introduce these three-dimensional (3D) heterogeneities in the surface formulation, models with several surfaces have been heuristically proposed [6].
This publication supplements two previous papers on two-layer atrial models.The biomedical relevance of the two-layer model was widely numerically investigated in [4] whereas a two-layer model of human atria was presented in [10].The present paper addresses the consistency of these models and rigorously proves, under simplifying assumptions on the fiber orientation near the boundary, the mathematical convergence of the mono-and two-layer models towards averaged 3D solutions when the tissue thickness goes to zero.
The paper is organized as follows.Section 2 states the problem.Section 3 formalizes basic preliminary results that are used throughout the paper.Section 4 first presents a formal derivation of a second order surface model similar (but one order more accurate) to the usual surface models.The error estimate that justifies the formal expansion is then rigorously proved.In section 5, the averages of this second order model on each layer are shown to solve a set of two coupled surface monodomain equations.The error estimate of this two-layer model, respectively, to the transmural averages of the 3D model, is rigorously proved to be of order 3 .Section 6 is devoted to the numerical illustrations of the convergence theorems.
2. Problem statement.We introduce a 3D slab of tissue Φ = φ × (−h, h), where φ is an open subset of R 2 and 2h is its thickness.We suppose that histological homogeneities allow us to split the tissue into two superimposed layers Φ (k) , such that Φ (1) = φ × (0, h) and Φ (2) = φ × (−h, 0).The monodomain model depicts the evolution of the transmembrane potential u (k) h (t, x, z) ∈ R of the cardiac tissue and m unknowns w (k) h (t, x, z) ∈ R m that describe its electrophysiological state, where t is the time, and (x, z) is the position in φ × (−h, h).The indices k reflect that these unknowns are related to one tissue layer.A model of cardiac electrophysiology, formalized by two functions (f, g) of R × R m , couples u (k) h and w (k) h .For k = 1, 2 and t > 0, the monodomain model reads

.2)
A and C represent the myocytes surface-to-volume ratio and membrane capacitance.The electrical diffusion tensor is decomposed into a two-dimensional (2D) tensor σ (k) (x) ∈ L ∞ (ω) and a transmural conductivity σ (k) 3 (x) on each layer k, that only depend on the longitudinal variable x and are homogeneous along the transmural variable z.We suppose that there exist two real number 0 < σ < σ such that σ ≤ σ (k) 3 ≤ σ and σ|χ| 2 ≤ σ (k) χ • χ ≤ σ|χ| 2 for all χ ∈ R 2 .For the sake of simplicity, we further assume that there exists an open subset φ φ such that σ (1) = σ (2)  on φ \ φ .In other words, we suppose that the fiber distribution is homogeneous through the whole tissue thickness near its boundary.This assumption does not contradict histological knowledge of atrial fibers orientation [8].A similar assumption was discussed in [19].
A dimensionless version of the system (2.1)-(2.2) is necessary for its asymptotic analysis when the tissue thickness vanishes.Such a dimensionless model, including a dimensional study, has been precisely introduced in [4].We note the dimensionless space Ω = Ω (1) ∪ Ω (2) with Ω (k) = ω × Γ (k) , Γ (1) = (0, 1), Γ (2) = (−1, 0), and ω the dimensionless version of φ.Still noting the dimensionless coordinates (t, x, z), the dimensionless problem reads, for k = 1, 2 and t > 0, The parameters α, β are dimensionless parameters that gather information on the balance of different physical characteristics, e.g., A, C, max x∈R,y∈R m f (t, x, y), etc.The aspect ratio between the thickness of the tissue and its length is supposed to be small.This system is supplemented by the boundary and transmission conditions 3 ∂ z u (2) , u (1) = u (2) in ω × {0}, (2.6) where n is the outward normal to ω on its boundary.To simplify the situation, we assume that the initial data are independent of k ∈ {1, 2} and of z ∈ (−1, 1), specifically, with u 0 (x) and w 0 (x) given functions of x ∈ ω: Remark 1.To relax the assumption of a homogeneous initial condition, we can consider transmurally nonhomogeneous initial conditions such as u (k) (0, x, z) = u 0 0 (x)+ 2 u k,0 1 (x, z), where the couple (u 0 0 , u k,0 1 ) satisfy the boundary and transmission conditions (2.5)- (2.6).We can, alternatively, add a boundary layer in time near t = 0 to reach a situation where the following results are valid.

Technical preliminaries.
Along with the proofs of Theorems 1, 2, and 3, we will need the following elementary lemmas.

Lemma 1. Consider a function
dt, the integral remainder in the Taylor expansion.Denoting by ∇ 2 h the Hessian matrix of h, we have Lemma 3 (Gronwall).Suppose that y(t) ≥ 0 is a C 1 function of t > 0, and h(t) ≥ 0, d 1 (t) ≥ 0 are functions in L 1 loc ((0, +∞)) and that d 2 (t) in L 2 loc ((0, +∞)) is such that, for all t > 0, there exists for all t > 0, and y(0) = 0, then Proof.Let us prove Lemmas 1 to 3. Lemmas 1 and 2 are straightforward derivations of Taylor expansions with integral remainders.
Lemma 3 is a direct consequence of the usual inequality of Gronwall y(t) ≤ exp(c 0 t)y(0) We obtain the first inequality because y(0) = 0 and t 0 exp(c 0 (t − s))ds ≤ 1 c0 exp(c 0 t) and exp(c 0 (t − s)) ≤ exp(c 0 t) for 0 ≤ s ≤ t.Hence, we integrate this inequality to obtain The final result is a direct consequence of the integrated inequality: 4. The asymptotic model with one layer.4.1.Formal derivation.We consider the following expansion of (u We introduce this expansion in the system (2.3)-(2.4)and identify the coefficients having the same order with respect to 2 .The coefficient of order 1/ 2 yields the equation ∂ zz u (k) 0 = 0 for k = 1, 2 together with the boundary condition (2.5), in the direction z, and the transmission conditions (2.6) for u (k) 0 .These equations easily show that u (2) 0 (t, x, z) := u 0 (t, x) is a function independent of z, defined for t ≥ 0 and a.e.x ∈ ω.
We then can write a second order approximation of the functions (f, g): We then get the following equation on u with the boundary and interface conditions (2.5) to (2.6) on u (k) 1 .Since all the terms of the ODE (4.2) are constant through Γ (1) ∪ Γ (2) (initial condition w 0 and source function g(u 0 (t, x), .)), the solutions w (k) 0 are independent of z and k for all times t ≥ 0 and have the same value that we denote by w 0 (t, x).
We denote by σ m = σ (1) +σ (2)   2 the arithmetic average of the conductivity matrices in both layers.We also introduce the notation 2 that will be used below.
Remark 3. Due to the definition of σ (k) on ω \ ω , b = 0 on that domain.This transmurally homogeneous distribution of fibers near the boundary, realistic from a physiological point of view [8,16], permits the derivation of the H 1 -type estimates in Theorems 1 and 2. A tangential distribution of the fibers, respectively, to the boundary, would also have been an acceptable simplification to derive this estimate.

The transmission condition on u (k)
1 on ω yields c (1) = c (2) := c so that through each layer k = 1, 2 and find that ū(k) As a consequence, we have the relations 3 σ (1) 3 where ū1 denotes the average of u through the whole thickness of the tissue and 3 +σ (2) 3 is the harmonic average of the transverse conductivities.Hence, c can be found if we know ū1 .We also denote by w1 the average of w through the whole thickness of the tissue, defined by w1 := 1 (•, z)dz.We now identify the coefficients of order 2 in the expansion of (u (k) , w (k) ).We first get the equations with the boundary and transmission conditions (2.5) to (2.6) on u (k) 2 .For k = 1, 2 and for each (x, z) ∈ Ω (k) , (4.10) is a first order linear Cauchy problem on w We note that the function A(t) is independent of z and the functions z → B (k) (t) (k = 1, 2) are polynomial functions for all t ≥ 0 and a.e.x ∈ ω.As a consequence, the functions z → w (k) 1 (t, x, z) are also polynomials for all t ≥ 0 and a.e.x ∈ ω.Remark 4. We additionally remark that w satisfies the transmission conditions (2.6).
Afterwards, we integrate again (4.9) for z ∈ Γ (k) , add the resulting equations, and use the conditions (2.5) and (2.6) on u (k) 2 .We remark that , and finally obtain the following equations for (ū 1 , w1 ): with the boundary and initial conditions

Definition of the first order solution and associated asymptotic problem.
We now redefine (u 0 , w 0 ) and (ū 1 , w1 ), which are no longer the coefficients of a formal asymptotic approximation, but the solutions to the bidimensional systems (4.4), (4.5) and (4.11), (4.12), respectively, for t > 0 and x ∈ ω, and with the boundary and initial conditions (4.6) and (4.13).In (4.11) and (4.12), the function b is defined on (0, +∞) × ω by b = div x (σ d ∇ x u 0 ).Afterwards, we define the 3D functions u 3 σ (2) 3 and the function w is the solution to the system of equations (4.10) for k = 1, 2 that also reads, with A = ∂ 2 g(u 0 , w 0 ) and Remark 5. Some straightforward computations show that, conversely, the averages in the thickness of u are exactly (ū 1 , w1 ), solutions of system (4.11)-(4.12).We also remark that (4.1) holds for u (k) 1 and u 0 .Finally, we check that the boundary conditions on u 0 and ū1 can be written k=1,2 σ (k) ∇ x u 0 • n = 0 and We note (ũ (k) , w(k) ), the first order approximate solution on (0, +∞) × Ω (k) : and we introduce (e (k) , f (k) ) as the corresponding errors with respect to the complete 3D solution (u

4.3.
Error estimate for the one-layer model.In Theorem 1 below, we prove that the errors e (k) and f (k) are bounded by 3 in L 2 (Ω (k) ) for all time t > 0, and that e (k) is also bounded by 3 in L 2 (0, t; H 1 (Ω (k) )) also for all time t > 0.
Theorem 1 (error estimates for the one-layer model).We assume that the functions f and g are C 2 (R × R m ) and that the solutions (u (k) , w (k) ), (u 0 , w 0 ), and uniformly in time and .Namely, we require that there exists M > 0 such that, for all t > 0 and and for all > 0, for all t > 0, and We now define the functions ψ (1) )dζ where the functions φ (k) are given by (4.21) below.We assume that . Then, we have the following estimates, for all 0 < ≤ 1, t > 0, and k = 1, 2: where f is defined below.Remark 6.Let us exhibit a class of source functions, which includes the Hodgkin-Huxley model and several of its adaptations to the cardiac cell, that ensures a uniform bound on (u (k) , w (k) ) with respect to , as required in Theorem 1.Let f and g be defined by , where h k , w k,∞ , and τ k are regular functions.In the electrophysiology literature, h k and τ k are strictly positive and w k,∞ belongs to (0, 1).Then it can be shown that the domain [min u k , max u k ] × [0, 1] m is invariant for the corresponding monodomain problem [17].Namely, for a given > 0, the diffusion tensor of problem Comparison arguments are then used to show that, for any initial condition belonging to for any t > 0. We can see that the invariant domain only depends on the source function parameters u k , which are independent of .The bounds on (u (k) , w (k) ) are then uniform in .In section 6, we use the Beeler-Reuter model which includes an intracellular concentration of calcium as state variable that does not fit the previous family.But it belongs to a class of ionic model whose boundedness has been studied for the bidomain model in [19].Again, the bounds exhibited in that reference only depend on the source function parameters, which are independent of > 0.
The functions φ (k) play the role of the functions 2 in the formal expansion of the previous section.The functions ψ (1) We average (4.21) in the whole thickness; use the result from Remark 5 and the tricks used to establish (4.11) to prove that for all t > 0 and a.e.x ∈ ω, because ū1 and w1 are solutions to the system (4.11),(4.12).It shows the transmission condition σ  k) , u (k) ), (u 0 , w 0 ), and (u 1 ) immediately gives the following equations: ).For k = 1, 2, we multiply (4.22) by e (k) and (4.23) by αf (k) , integrate on Ω (k) , and add the resulting equations in order to obtain the following energy estimate: where ) is the square of the total L 2 norm of the error and, with the notation To get the term A, we integrated by parts in the z direction with the boundary and transmission conditions (2.5)-(2.6)for e (k) , which nullified the trace on ω × {±1} and balanced the trace on ω × {0} of e (k) ψ (k) for k = 1, 2. We are going to control independently each term of the right-hand side of that estimate.Estimate on A. Using Young's inequality, we can see that

Estimate on B.
We remark that simple computations on the boundary conditions, using σ (k) ∇u (k) • n = 0 on ∂ω, give Each term of the right-hand side can be bounded independently on ∂ω.

Estimate on C. We define
and use Lemma 1 to check that, for k = 1, 2, (4.25) where Λ 0 and Λ 1 are defined in the statements of Theorem 1 and I f is defined in Lemma 1.We have used the fact that |(ũ for ≤ 1, because of the assumptions on (u 0 , w 0 ) and (u ).Now, we can write the last estimate: Gathering the estimates on A, B, and C, passing the z-derivative term of the A estimates to the left-hand side, and assuming that ≤ 1, i.e., 8 ≤ 6 , we get where the values c 0 , d 0 , and d 1 (t) are defined in the statement of the theorem.Lemma 3 then ends the proof.
We then focus on error estimates on the z derivative, which will be used later to prove the convergence of the two-layer model.
Theorem 2 (error estimates on the z derivative).Under the same assumptions as in Theorem 1 and assuming that k=1,2 ∂ z ψ (k) (s) 2 L 2 (Ω (k) ) belongs to L 1 loc (R + ), we have the following estimates for all 0 < ≤ 1, t > 0, and k = 1, 2, , and where c 2 = 2(Λ 0 +Λ 1 ), d 2 = 2ασM 4 Λ 1 |ω|, and integrate on Ω (k) , and add the resulting equations in order to obtain where y z (t) = α k=1,2 σ and To get this estimate, we used two ingredients.First, we integrated by parts on the z direction with the boundary conditions (2.5) and the transmission conditions (2.6) for e (k) and f (k) , which nullified the trace on ω × {±1} and balanced the trace on ω × {0} Second, we computed successive integration by parts along z and x, with transposition of the differential operators.We again control independently each term of the right-hand side of (4.31).
Estimate on A z .A Young's inequality first shows that (4.32) Estimate on B z .Recalling that σ We integrated twice by parts in the z direction with the boundary and transmission conditions (2.5)-(2.6)for e (k) , we used the boundary condition on u (k) on ∂ω, and the fact that u 0 is constant along z.Recalling (Remark 3) that b = 0 on ω \ ω , we know that σ (k) ∇b • n = 0 a.e. on ∂ω.We then finally have B z = 0. Estimate on C z .We first remark that and that (∂ z u We complete the estimate with (4.33) where we used some Young's inequalities, the definition of y z (t), and the inequality Using estimates (4.32) and (4.33), the final estimate then reads where the constants c 2 , d 2 , and the function d 3 (t) are given in the statement of Theorem 1.We conclude with Lemma 3.
Remark 7. The final error estimate is of order 3 , whereas a continuation of the formal study in the beginning of the section 4 would show that the coefficient of order 3 vanishes, and that the convergence order is of order 4 .This is a usual fact when estimating the error of the asymptotic expansion by a L 2 energy method; see [11].
5. The asymptotic model with two layers.
5.1.Equations on the averages through the thickness of each layer.In this section, we prove that the averages in z ∈ Γ (k) for each of the two layers of the solution (u (k) , w (k) ) verifies a system of two surface monodomain equations, linearly coupled, up to an error of order 3 .It is the basis of our two-layer model.We recall, for all integrable function z → h(z) on Γ (k) , the notation h = Γ (k) h(s)ds.
We integrate (4.16) for z ∈ Γ (k) and use the properties of u In particular, we can see that We derive (4.16) in z and obtain ∂ z e (k) .Now, according to the flux continuity of the transmission condition (2.6), we can define the residual Su (t, x) := σ 3 ∂ z u (2) (t, x, 0) = 2 b(t, x)+Su (t, x).We are ready to integrate (2.3) and (2.4) for z ∈ Γ (k) : As in the proof of Theorem 1, we denote by ) the nonlinear errors.Afterwards, we substitute b according to (5.1), so that the averages (ū (k) , w(k) ) are solutions to the system of equations for k = 1, 2 and setting

Definition of the two-layer model.
Consider some initial data û(k) 0 (x) and ŵ(k) 0 (x) defined on ω.The two-layer model is given by the following coupled system of 2D monodomain equations for k = 1, 2: It is a perturbation of the problem (5.2)-(5.3).We are going to show that this perturbation is actually small.

Error estimate for the two-layer model. We introduce the notations ê
Theorem 3 (error estimates for the two-layer model).Assume that, for k = 1, 2, we have û(k) 0 (x) = u 0 (x) and ŵ(k) 0 (x) = w 0 (x), and that, with the bound Under the assumptions of Theorem 1, we have the following estimates: for all 0 < ≤ 1, for all t > 0, where the functions k 4 and . We use the constants c 4 = 1 + 2Λ 0 and d 4 = 8αΛ 2  1 M 4 |ω|, and the functions d 5 (t) = 16Λ 2 0 αk 0 (t) 2 exp(c 0 t) and Proof.We closely follow the method presented in the previous proofs.We subtract (5.4) and (5.5) from equations (5.2) and (5.3) and obtain, for k = 1, 2, We multiply the equations by ê(k) and α f (k) , integrate on ω, sum for k = 1, 2, and finally obtain the energy estimate, The nonlinear terms read ).We note that, unlike the previous proofs, there is no trace residual on the boundary of ω, because of the boundary conditions on ū(k) and û(k) .Estimate on Â.With Young's inequality, we first remark that We then focus on estimating the first term: Now recall the definition of Su and note that, for all t > 0, for a.e.x ∈ ω, we have 3 ∂ zz e (1) (t, x, z)dz so that, because ∂ zz e (1) dz ∂ zz e (1) 2 dx.
Here, the choice of e (1) is arbitrary and could be replaced by e (2) .As a consequence, we also have ∂ z e (k) 2 dx with the transmission condition and Lemma 2. We finally have Last, recalling that ( ū(k) , w(k) ) = (u 0 , w 0 ) + 2 (ū For these three terms, the Lipschitz continuity of f indicates that e (k) , f (k) dz.
The final energy estimate reads (for 0 < ≤ 1) where )) 1/2 .The constants c 4 , d 4 , and the functions d 5 (t) are given in the statement of Theorem 3. Lemma 3 proves the result, because, for all t > 0 and 0 < < 1, according to Theorems 1 and 2, t 0 d 6 (s) 2 ds < C d6 (t) 6 .Remark 8.This theorem guarantees that the solution of the two-layer model converges toward the average in z by layer of the three-dimensional potential.Furthermore, the accuracy of the two-layer model is limited by the precision of the approximation of the transverse diffusion, i.e., of the asymptotic expansion of u (k) .
6.1.The test cases.We want to simulate the propagation of the activation front on a 3D slab of tissue Ω = ω × (−h, h) for various values of the thickness parameter h > 0.We take ω = (0, x 0 ) × (0, x 0 ) and define the conductivity matrices σ (k) in the layers Ω (k) such that the unitary eigenvector associated with the principal eigenvalue of each σ (k) is normal.(The fiber directions of both layers are perpendicular.)We have the choice between three models: 3D: The complete 3D equations (2.1) and (2.2) on the layers Ω (k) , which solutions are denoted by u h .1×2D: The order 0 term of the asymptotic expansion, u 0 and w 0 -discarding the corrector terms 2 (u 1 )-that solve the usual surface monodomain equations (4.4) and (4.5) on ω (independent on h).2×2D: The solution of the two-layer model, û h and ŵ(k) h , that solves (5.9) and (5.10) on ω with the coupling parameter σ 3 +σ (2) 3 .
The functions f and g are defined by the Beeler-Reuter ionic model [1], for it is computationally simple but retains the essential features of the cardiac action potential.The remaining parameters of the equations are given in Table 1.We will explore values of the thickness parameter h ranging from the typical diameter of a cardiomyocyte (tens of µm) up to large atrial thickness (less than a cm) as found in human tissues.An action potential is initiated by applying a transmembrane voltage of 20 mV during a period of 2 ms in the domain S = (0.49x 0 , 0.51x 0 ) × (0.49x 0 , 0.51x 0 ) for the surface models and in the cylinder S × (−h, h) for the 3D model.

Meshes, discretization, and resolution.
We consider a Cartesian grid of Ω with space steps Δx × Δx × Δz in the directions x, y, and z; and the subgrid with space steps Δx × Δx of ω.These Cartesian grids are split into a simplicial mesh of Ω and ω (that is, respectively, with tetrahedra and triangles), by splitting each square into 4 triangles and each hexahedron into 24 tetrahedra.In order to guarantee the convergence of the discrete solution, we take Δx = x 0 /200 = 5.e − 3 cm and Δz is chosen as in Table 2.This choice guarantees that there are at least 10 elements through the whole thickness of Ω (that is 2h) and Δz ≤ 4Δx.
The different equations are discretized by using the standard P1-Lagrange finite element method in space with diagonal mass lumping (due to the numerical quadrature rule) and the Rush-Larsen time-stepping scheme [15] with the fixed time step Δt = 0.05 ms.During the Rush-Larsen iteration, the diffusion terms are solved implicitly; so are the coupling terms in the coupled equations of the two-layer model.All the linear systems are solved with a conjugate gradient method, a Jacobi preconditioner, and a fixed tolerance equal to 1.e − 10. 2x2D (e) û( 1) In the asymptotic regime, all solutions are qualitatively very similar (Figures 1(a), 1(b), 1(e), 1(f) and 1(i)), as expected.In the nonasymptotic regime, the usual surface model clearly misses the important 3D interplay of diffusion and reaction between the two layers of tissue [4].This interplay induces the diamond-shaped wavefront seen on Figures 1(c)-1(d), which remains in our 2×2D model (Figures 1(g)-1(h)), but not in the usual 1×2D model (Figures 1(i)).Note that such a wavefront was observed experimentally in [20].Our 2×2D model clearly accounts for such phenomena.

Errors.
In order to quantify the distances from the 1×2D and 2×2D models, respectively, to the 3D one, and, accordingly, with the results from Theorems 1 and 3, we define the following errors E 0 and E 1 :  We expect that E 0 (h) = O(h 2 ) and E 1 (h) = O(h 3 ), although we presume that E 1 (h) = O(h 4 ) is the correct convergence rate (Remark 8).The error graphs (h, E 0 (h)) and (h, E 1 (h)) are displayed in Figure 2. As expected, we observe the convergence of both models, and a higher convergence rate for the two-layer model (E 1 (h)) than for the surface model (E 0 (h)), but only for h < 0.1 cm.Remark 9. We point out that the convergence graph in Figure 2 differs from the graph presented in [4], due to the different metrics that were used (we used L ∞ (0, T ; L 2 (ω)) relative errors in Figure 2 versus L 2 (ω) relative errors on activation maps in [4]).
Finally, we emphasize the very strong improvement of the computational cost for the 2×2D (25 to 70 times faster), and the usual 1×2D (70 to 170 times faster) models, respectively, to the 3D one.

Interest of the two-layer model.
A similar two-layer model was heuristically used in [6], but the rigorous mathematical derivation, the error estimates proved in Theorem 3, and the numerical study from section 6 are new.
The main advantage of the 2×2D model over the usual 1×2D model resides in the distinction of individual reaction terms 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 [20].The usual surface model, even enhanced by addition of the second order term 2 u (k) 1 would not easily account for these phenomena.Yet, we point out that these phenomena (dissociation, transmural gradients, complex propagation) are expected to play a major role in the initiation of arrhythmias.Hence, the two-layer model is a very efficient tool to investigate in depth the arrythmogenic mechanisms in the atria, providing a huge reduction of the computational load comparatively to the 3D model.

Limitations and perspectives.
Although being a strong qualitative enhancement of the usual 1×2D model, the 2×2D model could still be improved, for instance, by incorporating higher order terms in the Taylor expansions.The deriva-tion of a more general bidomain version of the two-layer model on more general geometries is also possible, following the work from [2].We also note that simplifying assumptions of this work can be relaxed; cf., e.g., Remark 1.
Our numerical study illustrates the results from sections 4 and 5 and we give errors in the L 2 norm.In [4], we investigate the 2×2D model by comparing its activation maps to the 3D ones for different geometries; this work addresses the 2×2D robustness to meaningful biomedical situations.Furthermore, we propose in [10] an anatomically and structurally accurate 2×2D model of human atria, including fiber and functional heterogeneities, which is able to address biomedical issues.We believe that it will become an appreciated tool to study the structure-to-function relations in the atria, from a medical point of view.
, x, z) = u where b = div x σ d ∇ x u 0 , and c = c(t, x) is an unknown function.We denote by ū

h
refer to the discrete solutions of the models at time t n = nΔt.

Table 1
Parameters of the equations.

Table 2
Meshes and numbers of degrees of freedom (#dof ). h