A computational bilayer surface model of human atria

A bilayer surface model of human atria is presented. This model allows inclusion of transmural heterogeneities while maintaining small computational costs associated with surface models. A physiological fibre structure is constructed on an imaged atrial geometry which includes the main structures of fibres in the right and left atria, including the transmural heterogeneities, and the main transseptal connexions. The validity of the corresponding bilayer mathematical model is assessed by comparing the case of two sheets of orthogonal fibres to a three dimensional slab of tissue. A comparison of bilayer and monolayer simulations of sinus propagation is finally proposed. This atrial bilayer model is a good tool for performing complex atrial simulations in pathological cases.


Introduction
Atrial computational models have proven their utility in fundamental or clinical research by giving insight into atrial fibrillation (AF) initiation and perpetuation, radio-frequency ablation, and the development of signal analysis techniques [1]. 1 As the electrical activity spreads across the atria, propagation patterns and P-wave duration are strongly influenced by atrial fibre structure [2,3]. Modeling the anisotropy of the atria is, therefore, a major step in the construction of atrial models. This construction can result from an imagebased process on ex-vivo hearts [3] or from a semi-automatic rule-based approach for acquired geometries of in-vivo atria [2,4].
Atrial three dimensional models are usually divided into two main categories : surface and volumetric models. Whereas the former propose a homogeneous transmural fibre distribution and computationally light simulations, the latter capture important transmural heterogeneities or layer dissociations as observed in animal models [5,6,7] in exchange for increased computational costs [1].
As a compromise between the surface two dimensional and volumetric three dimensional models, we propose a two layer computational model of the atria that can reproduce complex three dimensional patterns while keeping surface computational load. Such a two layer approach has also been used in [8] on idealized geometries, but has not been demonstrated to reproduce proper behavior.
We first describe the construction of the model. We then give some mathematical arguments to the validity of the bilayer approach. We finally compare simulations of a normal action potential (AP) propagation in a bilayer and a monolayer model.

Imaging the Atria
Bi-atrial geometry was acquired in vivo using ECG-gated contrast-enhanced multi-detector computed tomography (MDCT) in a patient with no history of atrial disorder. Imaging was performed using a 64-slice CT scanner (SO-MATOM Definition, Siemens Medical Solutions, Forchheim, Germany) during intravenous injection of a 120 mL bolus of iodinated contrast agent at the rate of 4 mL s −1 . ECG-gating was set for the acquisition window to occur at end-systole, when atrial motion is minimal. Acquisition parameters were: slice thickness 0.5 mm, tube Voltage 120 kV, maximum tube current 850 mAs, and gantry rotation time 330 ms. An imaging volume comprising the whole atria was reconstructed with a voxel size of 0.5x0.4x0.4 mm. Transaxial images were imported in DICOM format in a local database of the software OsiriX 3.6.1 (OsiriX foundation, Geneva, Switzerland). The endocardial surface was segmented automatically using region growth segmentation.
Volume rendering reconstructions were used to remove non atrial structures from the segmentation. Pulmonary veins were cut several centimetres from their ostia because myocardium is known to extend to the first centimeter beyond the pulmonary veins ostia. The resulting segmentation was used to compute a 3D mesh in vtk format using the software CardioviZ3D (IN-RIA Asclepios, Sophia Antipolis, France). The mesh was smoothed with weighted Laplacian smoothing and optimized with the software freeYams to obtain a mesh suitable for calculation that contained 361483 nodes and 715034 elements with a mean diameter of 0.492 mm.

Algorithms.
Two methods were used for fibres definition and mesh completion.
Fibre definition. The fibre structure was defined with a rule-based semiautomatic method [4] that was founded on two atrial characteristics : (1) atrial geometry can be divided in basic cylinder-like structures with circumferential fibre orientation. (2) Fibre arrangement is continuous. The atrial mesh was first manually decomposed in cylinder-like elements. The local circumferential structure was semi-automatically obtained by determining the cylindrical axis as an eigen vector of an inertia tensor. Finally, an inpainting algorithm extended continuously the fibre direction on the whole geometry.
Manual mesh completion. Small structures that were partially acquired, such as transseptal connexions, were manually constructed with the following workflow : the points of insertion of the new structure in the mesh were identified. Plane surfaces connecting those points were defined and meshed with the software Gmsh. The structure was then smoothed by a Laplacian method with the software MeshLab (3D-CoForm Project, Italy).

Main modeling options.
In order to describe the atria as a bilayer structure, we need to define some physiological a-priori modeling options, specially focusing on transmural heterogeneities.
Right Atrium (RA). According to histological descriptions of the atria [10], the RA can be seen as a monolayer structure, except in the appendage (RAA) and in the roof where the crista terminalis (CT) and the pectinate muscles (PM) bring thickness and fibres inhomogeneities. A first layer modeled the epicardial part of the RA including its main fibre structures : the circular fibres that surrounded the superior veina cava (SVC) crossed the CT and covered the RAA to blend with the circular bundle of the vestibule. A second layer was applied on the location of PM and the CT to include the longitudinal fibre direction along those structures in the endocardial part of the RA. A poor conduction zone was introduced along the left side of the sino-atrial node (SAN) as assessed in [12].
Left Atrium (LA). The left atrium (LA) fibres structure can be described as a two layer geometry. The main transmural heterogeneities are concentrated in the anterior wall and on the pulmonary veins (PV). The circumferential and oblique fibres of the anterior epicardium blend with the Bachmann's bundle whereas the endocardial fibres of the septoatrial bundle run from the vestibule to the dome of the LA. The fibre structure of the PV can show strong heterogeneities with non uniform arrangement and transmural discontinuities of fibres directions [11].
Transseptal connexions. The chambers of this atrial model were electrically insulated and connected by the three main interatrial connexions : the BB, the fossa ovalis (FO) and the coronary sinus (CS). The BB was segmented manually as a three-dimensional structure and then reduced to two surfaces : the first one corresponded to the posterior face of the BB and inserted into the septal wall of LA while the second one was constructed as a mean surface of the volumetric BB and blended into the septopulmonary bundle of the LA. The FO and the CS were identifiable after the segmentation step and connected manually.

Modeling the AP Propagation
The monodomain equation was used to describe the propagation of electrical conductivity where V m is the transmembrane voltage, σ i is the conducitivity tensor, β is the ratio of membrane surface to volume in which it is contained, C m is the membrane capacitance (1F cm −2 ), and I ion is the transmembrane current which is a function of voltage, time, and ζ, a set of state variables.
Mathematical assessment of the bilayer model validity. A bilayer expression of this mathematical model has been asymptotically derived in the limit of vanishing thickness in [9] from a bilayer volumetric monodomain model. That reads: . σ ′(k) and σ (k) 3 are respectively the diffusion coefficients in the layer tangent plane and in the transmural direction. h is the thickness of the layers.
Three dimensional, bilayer and monolayer finite element P 1 -Lagrange simulations have been performed on, respectively, a 3D slab and its 2D equivalent for various h, with the Beeler-Reuter model and a semi-implicit Euler scheme. The activation maps of the mean potential in the thickness of each simulation have been compared. Variations of the coupling coefficient σ c were uniformly applied to check the robustness of the model under physiological thickness, far from the asymptotic regime.
Bilayer atrial model. The ionic model was described by the Courtemanche human atrial action potential model [13]. The sodium conductance was increased by a factor of four from the published value to achieve realistic propagation speeds. σ i and β were set to (0.3, 0.04) S m −1 and 0.14 m −2 . The structural block zone close to the SAN [12] was implemented with a passive ionic model with a resting level -80 mV and an intracellular conductivity reduced by a factor of 100. The SAN is located between the block zone and the crista terminalis. Accordingly, a small region was stimulated there to simulate the site of first excitation.
The equations were solved using the finite element method. Each surface was composed of linear triangular elements while connections between surfaces were implemented as one dimensional finite elements. 134510 line elements were added to couple the surfaces together. The system was solved using the CARP cardiac simulator [14] with a fixed time step of 25 s using a Crank-Nicolson approach to solve eq. 1. On a desktop workstation with Intel(R) Xeon X5650 CPUs running at 2.67GHz, it took seven seconds to simulate one ms on ten cores, while on a desktop with an NVidia Quadro 4000 (256 cores) it took 18 seconds per millisecond running it on the GPU.

Fibres architecture
The main fibre structures as described in [10] are correctly reproduced by our model. In fig. 2, we can see that the differences between the endo and epicardial fibres orientations are mainly localized in the CT and PM in the RA, and in the anterior and posterior wall, together with the PVs in the LA, where the fibres in both layers can be orthogonal.

Simulations
Validity of the bilayer approach. The bilayer mathematical model has been shown to have second-order convergence towards the volumetric model when h −→ 0. For physiological h, a correction can be brought to the coupling value σ c so that the accuracy of the bilayer surface model is preserved : complex depolarisation patterns can be obtained with a relative error of 3% respectively to 3D simulations ( figure 3). We refer to [9] for a deeper study of this model.
Bi-atrial simulations. The propagation of electrical activity showed a normal pathway: it first propagated from SAN to the BB along the CT.

Discussion and conclusions
Clinical uses of atrial modeling or acurate atrial models calibrated for fundamental research need personalized fibres orientations, and then efficient methods of construction. The work-flow presented here is still time-consuming and has to be improved in view of such a systematic use. Other methods that are being developped for 3D models [2] could be adapted to bilayer surface models. Statistical methods of registration to transport this model to personalized geometries could be another way to shorten this process [15]. Modeling assumptions were set prior to model construction. The bilayer representation of some zones (LA,CT,PM) and the monolayer representation of others (RA) were guided by histological descriptions. If more complex transmural heterogeneities are observed, the model can be locally improved by adding layers or by simulating small entities with a volumetric model.
A normal propagation was obtained. The main qualitative characteristics of an atrial normal AP propagation were correctly reproduced : the transseptal pathway and the endpoint of the sinus wave are in a physiological range. The total activation time is larger than those reported by the literature [16] but this can be tuned by further adjusting a combination of sodium channel conductances and tissue conductivities.
The bilayer model demonstrated its capacity to capture some perturbations of action potential propagation induced by transmural heterogeneities; that represents an important improvement over monolayer surface models. We could show that those perturbations were reliable on idealized geometries, indicating that this model could be efficiently used on physiological geometries. The demonstration of the potential of the bilayer approach could be completed with a comparison between a volumetric model and its equivalent bilayer surface model. This bi-atrial bilayer model proved to be a good trade-off between surface and volumetric models; it reproduced transmural heterogeneities that were specific to volumetric models with a small computational load. It is an